Dante
Dante

Reputation: 31

Intermediate points between 2 geographic coordinates

I am trying to develop an algorithm that involves normalizing GPS coordinates (latitude/longitude). That means, that being given two points A (lat1,lon1) and B(lat2,lon2) I would like to insert a point C that is linear with AB (same arc) and is placed at a specific distance from A and B (eg: A to B distance is 0.5km and I want point C to be at 0.1 km from A, on the AB arc). How can I calculate the coordinates for point C? For the purpose given, it is enough to approximate Earth as a perfect spherical object. I have found this article, but it gives the formula for midpoint only (and I don't fully understand it, in order to adapt). midpoint between two latitude and longitude Thank you.

Edit: I tried this but it gives wrong answers

public static void normalizedPoint(double lat1, double lon1, double lat2, double lon2, double dist){
        double constant=Math.PI/180;
        double angular = dist/6371;
        double a = Math.Sin( 0* angular )/Math.Sin(angular);
        double b = Math.Sin(1*angular)/Math.Sin(angular);
        double x = a * Math.Cos(lat1) * Math.Cos(lon1) + b * Math.Cos(lat2) * Math.Cos(lon2);
        double y = a * Math.Cos(lat1) * Math.Sin(lon1) + b * Math.Cos(lat2) * Math.Sin(lon2);
        double z = a * Math.Sin(lat1) + b * Math.Sin (lon2);
        double lat3 = Math.Atan2(z, Math.Sqrt( x*x + y*y ));
        double lon3 = Math.Atan2(y, x);
        Console.WriteLine(lat3/constant + " " + lon3/constant );
    }

As far as I understood the original formulas this should return one of the 2 original points, but it does not(because the fraction used is 1). Also the variable dist is the distance from the 2 points and is properly calculated (checked with the same website).

Edit 2: I am providing as inputs coordinates for 2 geographic points (lat1, lon1, lat2 lon2) and the distance between them. I'm trying to get an intermediary point (lat3,lon3).

Upvotes: 3

Views: 3321

Answers (3)

Tasdiqul Islam
Tasdiqul Islam

Reputation: 61

def get_intermediate_point(lat1 , lon1 , lat2 , lon2 , d):
    constant = np.pi / 180
    R = 6371
    φ1 = lat1 * constant
    λ1 = lon1 * constant
    φ2 = lat2 * constant
    λ2 = lon2 * constant
    y = np.sin(λ2-λ1) * np.cos(φ2);
    x = np.cos(φ1)*np.sin(φ2) -  np.sin(φ1)*np.cos(φ2)*np.cos(λ2-λ1)
    θ = np.arctan2(y, x)
    brng = (θ*180/np.pi + 360) % 360;  #in degrees
    brng = brng * constant

    φ3 = np.arcsin( np.sin(φ1)*np.cos(d/R ) + np.cos(φ1)*np.sin(d/R )*np.cos(brng) )
    λ3 = λ1 + np.arctan2(np.sin(brng)*np.sin(d/R )*np.cos(φ1),  np.cos(d/R )-np.sin(φ1)*np.sin(φ2));

    return φ3/constant , λ3/constant

Upvotes: 0

Witcher
Witcher

Reputation: 1140

Not sure if the original author found some answer, but since I had similar problem and developed working solution, I think it would be good to post it here.

The Problem

Having two geographical points, A and B, find intermediate point C which lies exactly on the direct way from A to B and is N kilometers far from A (where N is less than distance between A and B, otherwise C = B).

My Context

I was developing small pet project based on microservices architecture. The idea was to launch missile from given deployment platform (point A) to chosen target location (point B). I had to create some kind of simulator that sends some messages about current missile Geo location after it is launched, so I had to find those intermediate points between A and B somehow.

Solution Context

Eventually, I developed C# based solution based on this great web page - https://www.movable-type.co.uk/scripts/latlong.html.

That web page has all the explanations, formulas and JavaScript code at the bottom. If you are not familiar with the C#, you can use their JavaScript implementation.

My C# Implementation

Your input is Location A, Location B and the distance.

  1. You need to find bearing from A to B (see 'Bearing' section on that site)
  2. You need to find position C from A having bearing (see 'Destination point given distance and bearing from start point' on that site)

The Code

I have working solution as part of my pet project and it can be found here - https://github.com/kakarotto67/mlmc/blob/master/src/Services/MGCC.Api/ChipSimulation/CoordinatesHelper.cs. (Since original class is subject to change in the future, you might need to refer to this gist - https://gist.github.com/kakarotto67/ef682bb5b3c8bd822c7f3cbce86ff372)

Usage

// 1. Find bearing between A and B    
var bearing = CoordinatesHelper.FindInitialBearing(pointA, pointB);

// 2. Find intermediate point C having bearing (above) and any distance in km
var pointC = CoordinatesHelper.GetIntermediateLocation(pointA, bearing, distance);

I hope somebody will find this helpful.

Upvotes: 3

Damien_The_Unbeliever
Damien_The_Unbeliever

Reputation: 239636

As I point out in an answer on the linked to question, you need to change all of your inputs to use radians rather than degrees.

I believe you also had an error for z where you used lon2 rather than lat2.

With those corrections, I get the answer you're seeking:

    public static void normalizedPoint(double lat1, double lon1,
                                       double lat2, double lon2,
                                       double dist)
    {
        double constant = Math.PI / 180;
        double angular = dist / 6371;
        double a = Math.Sin(0 * angular) / Math.Sin(angular);
        double b = Math.Sin(1 * angular) / Math.Sin(angular);
        double x = a * Math.Cos(lat1* constant) * Math.Cos(lon1* constant) + 
                   b * Math.Cos(lat2* constant) * Math.Cos(lon2* constant);
        double y = a * Math.Cos(lat1* constant) * Math.Sin(lon1* constant) + 
                   b * Math.Cos(lat2* constant) * Math.Sin(lon2* constant);
        double z = a * Math.Sin(lat1* constant) + b * Math.Sin(lat2* constant);
        double lat3 = Math.Atan2(z, Math.Sqrt(x * x + y * y));
        double lon3 = Math.Atan2(y, x);
        Console.WriteLine(lat3 / constant + " " + lon3 / constant);
    }

Of course, the above can be vastly simplified by only converting angles ones, avoiding repeated calculations of the same Sin/Cos values, etc.

Calling:

normalizedPoint(47.20761, 27.02185, 47.20754, 27.02177, 1);

I get the output:

47.20754 27.02177

Upvotes: 2

Related Questions