7

I am currently using php to try and write code that will take a start point and end point latitude and longitude coordinate pair in decimal degrees and calculate the latitude and longitude of every vertex at a given interval in meters along the path between the two points.

For example:

I have start coordinates of latitude 43.97076 and longitude 12.72543 and I have an endpoint of latitude 43.969730 and longitude 12.728291. Using mathematical calculations (e.g. haversine formula) written in php, I would retrieve the coordinates between those two points at every 1 meter interval starting at the first coordinate and ending at the end coordinate. I welcome a solution written in any language like python, javascript, etc for an example.

  • Could you provide a little more information like what software or programming language you are using? – PMK Aug 10 '15 at 14:19
  • php, but answers in other programming languages are welcome too. – nulll Aug 10 '15 at 14:23
  • You say, precision in meters. Does this mean you want the starting point and a point at every meter along the line until the stop point? – JasonT Aug 10 '15 at 16:05
  • Yes! Exactly... – nulll Aug 10 '15 at 16:16
  • 2
    You will get better replies if you show your attempts to solve the problem (i.e. show your code) and are more specific about the issues. – MappaGnosis Aug 10 '15 at 16:45
9

This is quick and dirty. Most of this was adapted from the following website:

http://www.movable-type.co.uk/scripts/latlong.html

I haven't tested this too much but it seemed to work after initial testing. It will return a list of lat, long coordinate pairs along a line at a specified interval. I wrote it in python since I don't know php very well.

from math import cos, sin, atan2, sqrt, radians, degrees, asin, modf

def getPathLength(lat1,lng1,lat2,lng2):
    '''calculates the distance between two lat, long coordinate pairs'''
    R = 6371000 # radius of earth in m
    lat1rads = radians(lat1)
    lat2rads = radians(lat2)
    deltaLat = radians((lat2-lat1))
    deltaLng = radians((lng2-lng1))
    a = sin(deltaLat/2) * sin(deltaLat/2) + cos(lat1rads) * 
        cos(lat2rads) * sin(deltaLng/2) * sin(deltaLng/2)
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    d = R * c

    return d

def getDestinationLatLong(lat,lng,azimuth,distance):
    '''returns the lat an long of destination point 
    given the start lat, long, aziuth, and distance'''
    R = 6378.1 #Radius of the Earth in km
    brng = radians(azimuth) #Bearing is degrees converted to radians.
    d = distance/1000 #Distance m converted to km

    lat1 = radians(lat) #Current dd lat point converted to radians
    lon1 = radians(lng) #Current dd long point converted to radians

    lat2 = asin( sin(lat1) * cos(d/R) + cos(lat1)* sin(d/R)* cos(brng))

    lon2 = lon1 + atan2(sin(brng) * sin(d/R)* cos(lat1), 
           cos(d/R)- sin(lat1)* sin(lat2))

    #convert back to degrees
    lat2 = degrees(lat2)
    lon2 = degrees(lon2)

    return[lat2, lon2]

def main(interval,azimuth,lat1,lng1,lat2,lng2):
    '''returns every coordinate pair inbetween two coordinate 
    pairs given the desired interval'''
    coords = []
    d = getPathLength(lat1,lng1,lat2,lng2)
    remainder, dist = modf((d / interval))
    counter = 1.0
    coords.append([lat1,lng1])
    for distance in xrange(1,int(dist)):
        c = getDestinationLatLong(lat1,lng1,azimuth,counter)
        counter += 1.0
        coords.append(c)
    counter +=1
    coords.append([lat2,lng2])
    return coords

if __name__ == "__main__":
    #point interval in meters
    interval = 1
    #direction of line in degrees
    azimuth = 279.91
    #start point
    lat1 = 41.97809
    lng1 = -91.65573
    #end point
    lat2 = 41.97834
    lng2 = -91.65767
    coords = main(interval,azimuth,lat1,lng1,lat2,lng2)
    print coords

I hope this helps. sorry if not. :)

Ok, I added the ability for it to calculate the azimuth for you. I tested with the same location you were using and it turned out perfectly.http://www.darrinward.com/lat-long/?id=705140

import math

def getPathLength(lat1,lng1,lat2,lng2):
    '''calculates the distance between two lat, long coordinate pairs'''
    R = 6371000 # radius of earth in m
    lat1rads = math.radians(lat1)
    lat2rads = math.radians(lat2)
    deltaLat = math.radians((lat2-lat1))
    deltaLng = math.radians((lng2-lng1))
    a = math.sin(deltaLat/2) * math.sin(deltaLat/2) + math.cos(lat1rads) * math.cos(lat2rads) * math.sin(deltaLng/2) * math.sin(deltaLng/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = R * c
    return d

def getDestinationLatLong(lat,lng,azimuth,distance):
    '''returns the lat an long of destination point 
    given the start lat, long, aziuth, and distance'''
    R = 6378.1 #Radius of the Earth in km
    brng = math.radians(azimuth) #Bearing is degrees converted to radians.
    d = distance/1000 #Distance m converted to km
    lat1 = math.radians(lat) #Current dd lat point converted to radians
    lon1 = math.radians(lng) #Current dd long point converted to radians
    lat2 = math.asin(math.sin(lat1) * math.cos(d/R) + math.cos(lat1)* math.sin(d/R)* math.cos(brng))
    lon2 = lon1 + math.atan2(math.sin(brng) * math.sin(d/R)* math.cos(lat1), math.cos(d/R)- math.sin(lat1)* math.sin(lat2))
    #convert back to degrees
    lat2 = math.degrees(lat2)
    lon2 = math.degrees(lon2)
    return[lat2, lon2]

def calculateBearing(lat1,lng1,lat2,lng2):
    '''calculates the azimuth in degrees from start point to end point'''
    startLat = math.radians(lat1)
    startLong = math.radians(lng1)
    endLat = math.radians(lat2)
    endLong = math.radians(lng2)
    dLong = endLong - startLong
    dPhi = math.log(math.tan(endLat/2.0+math.pi/4.0)/math.tan(startLat/2.0+math.pi/4.0))
    if abs(dLong) > math.pi:
         if dLong > 0.0:
             dLong = -(2.0 * math.pi - dLong)
         else:
             dLong = (2.0 * math.pi + dLong)
    bearing = (math.degrees(math.atan2(dLong, dPhi)) + 360.0) % 360.0;
    return bearing

def main(interval,azimuth,lat1,lng1,lat2,lng2):
    '''returns every coordinate pair inbetween two coordinate 
    pairs given the desired interval'''

    d = getPathLength(lat1,lng1,lat2,lng2)
    remainder, dist = math.modf((d / interval))
    counter = float(interval)
    coords = []
    coords.append([lat1,lng1])
    for distance in xrange(0,int(dist)):
        coord = getDestinationLatLong(lat1,lng1,azimuth,counter)
        counter = counter + float(interval)
        coords.append(coord)
    coords.append([lat2,lng2])
    return coords

if __name__ == "__main__":
    #point interval in meters
    interval = 1.0
    #direction of line in degrees
    #start point
    lat1 = 43.97076
    lng1 = 12.72543
    #end point
    lat2 = 43.969730
    lng2 = 12.728294
    azimuth = calculateBearing(lat1,lng1,lat2,lng2)
    print azimuth
    coords = main(interval,azimuth,lat1,lng1,lat2,lng2)
    print coords

I think I have fixed the bugs from the original and it should work (including change the interval). Let me know how it turns out.

  • I succesfuly ported your code to php, obtaining the same results, but if i display points on a map, it seems that coordinates are going in the opposite direction, Look at the example darrinward.com/lat-long/?id=704976 – nulll Aug 10 '15 at 17:57
  • I forgot, You need to reverse the direction azimuth. – PMK Aug 10 '15 at 18:01
  • Updated and working now. – PMK Aug 10 '15 at 20:09
5

I did port the Python code to php:

<?php

function parseFloat($ptString) { 
            if (strlen($ptString) == 0) { 
                    return false; 
            } 

            $pString = str_replace(" ", "", $ptString); 

            if (substr_count($pString, ",") > 1) 
                $pString = str_replace(",", "", $pString); 

            if (substr_count($pString, ".") > 1) 
                $pString = str_replace(".", "", $pString); 

            $pregResult = array(); 

            $commaset = strpos($pString,','); 
            if ($commaset === false) {$commaset = -1;} 

            $pointset = strpos($pString,'.'); 
            if ($pointset === false) {$pointset = -1;} 

            $pregResultA = array(); 
            $pregResultB = array(); 

            if ($pointset < $commaset) { 
                preg_match('#(([-]?[0-9]+(\.[0-9])?)+(,[0-9]+)?)#', $pString, $pregResultA); 
            } 
            preg_match('#(([-]?[0-9]+(,[0-9])?)+(\.[0-9]+)?)#', $pString, $pregResultB); 
            if ((isset($pregResultA[0]) && (!isset($pregResultB[0]) 
                    || strstr($preResultA[0],$pregResultB[0]) == 0 
                    || !$pointset))) { 
                $numberString = $pregResultA[0]; 
                $numberString = str_replace('.','',$numberString); 
                $numberString = str_replace(',','.',$numberString); 
            } 
            elseif (isset($pregResultB[0]) && (!isset($pregResultA[0]) 
                    || strstr($pregResultB[0],$preResultA[0]) == 0 
                    || !$commaset)) { 
                $numberString = $pregResultB[0]; 
                $numberString = str_replace(',','',$numberString); 
            } 
            else { 
                return false; 
            } 
            $result = (float)$numberString; 
            return $result; 
}   

function xrange($start, $limit, $step = 1) {
    if ($start < $limit) {
        if ($step <= 0) {
            throw new LogicException('Step must be +ve');
        }

        for ($i = $start; $i <= $limit; $i += $step) {
            yield $i;
        }
    } else {
        if ($step >= 0) {
            throw new LogicException('Step must be -ve');
        }

        for ($i = $start; $i >= $limit; $i += $step) {
            yield $i;
        }
    }
}

function modf($x) {
    $m = fmod($x, 1);
    return [$m, $x - $m];
}


function getPathLength($lat1,$lng1,$lat2,$lng2)
{
    $R = '6371000'; //# radius of earth in m
    $lat1rads = deg2rad($lat1);
    $lat2rads = deg2rad($lat2);
    $deltaLat = deg2rad(($lat2 - $lat1));
    $deltaLng = deg2rad(($lng2 - $lng1));
    $a = sin($deltaLat/2) * sin($deltaLat/2) + cos($lat1rads) * cos($lat2rads) * sin($deltaLng/2) * sin($deltaLng/2);
    $c = 2 * atan2(sqrt($a), sqrt(1 - $a));
    $d = $R * $c;
    return $d;
}

function getDestinationLatLong($lat,$lng,$azimuth,$distance){

    $R = 6378.1; //#Radius of the Earth in km
    $brng = deg2rad($azimuth); #Bearing is degrees converted to radians.
    $d = $distance / 1000; #Distance m converted to km
    $lat1 = deg2rad($lat); #Current dd lat point converted to radians
    $lon1 = deg2rad($lng); #Current dd long point converted to radians
    $lat2 = asin(sin($lat1) * cos($d/$R) + cos($lat1)* sin($d/$R)* cos($brng));
    $lon2 = $lon1 + atan2(sin($brng) * sin($d/$R)* cos($lat1), cos($d/$R)- sin($lat1)* sin($lat2));
    #convert back to degrees
    $lat2 = rad2deg($lat2);
    $lon2 = rad2deg($lon2);

    return [$lat2, $lon2];  
}

function calculateBearing($lat1,$lng1,$lat2,$lng2){
   // '''calculates the azimuth in degrees from start point to end point'''
    $startLat = deg2rad($lat1);
    $startLong = deg2rad($lng1);
    $endLat = deg2rad($lat2);
    $endLong = deg2rad($lng2);
    $dLong = $endLong - $startLong;
    $dPhi = log(tan($endLat / 2 + pi() / 4) / tan($startLat / 2 + pi() / 4));
    if (abs($dLong) > pi()) {
        if ($dLong > 0) {
            $dLong = -(2 * pi() - $dLong);
        } else {
            $dLong = 2 * pi() + $dLong;
        }
    }
    $bearing = (rad2deg(atan2($dLong, $dPhi)) + 360) % 360;
    return $bearing;
}

function main($interval, $azimuth, $lat1, $lng1, $lat2, $lng2) 
{
    $d = getPathLength($lat1, $lng1, $lat2, $lng2);
    $rapydUnpack = modf($d / $interval);
    $remainder = $rapydUnpack[0];
    $dist = $rapydUnpack[1];
    $counter = parseFloat($interval);
    $coords = [];
    array_push($coords, [ $lat1, $lng1 ]);

    $xRange = xrange(0, intval($dist));

    foreach ($xRange as $rapydIndex => $value) 
    {
    //print $value;
        $distance =$value;
        $coord = getDestinationLatLong($lat1, $lng1, $azimuth, $counter);
        $counter = $counter + parseFloat($interval);
        array_push($coords, $coord);
    }
    array_push($coords, [ $lat2, $lng2]);

    return $coords;
}

  #point interval in meters
    $interval = 10;
    #direction of line in degrees
    #start point
    $lat1 = 43.97076;
    $lng1 = 12.72543;
    #end point
    $lat2 = 43.969730;
    $lng2 = 12.728294;
    $azimuth = calculateBearing($lat1,$lng1,$lat2,$lng2);
    //print $azimuth;
    $coords = main($interval,$azimuth,$lat1,$lng1,$lat2,$lng2);
    print_r($coords);

?>
  • This code seems to just hang – Barry Chapman Dec 7 '18 at 15:06
1

Thanks, @Peter for the code!
I converted it to NodeJS

/**
 * All lat-long specific operation define here
 *
 * http://www.movable-type.co.uk/scripts/latlong.html
 * https://gis.stackexchange.com/questions/157693/getting-all-vertex-lat-long-coordinates-every-1-meter-between-two-known-points
 *
 */
var _ = required('lodash');

function LatLong() {}

/**
 * Calculates the distance between two lat, long coordinate pairs
 * @param lat1
 * @param lng1
 * @param lat2
 * @param lng2
 *
 * @return integer 
 */
LatLong.prototype.getPathLength = function(lat1, lng1,lat2, lng2)
{
    var lat1rads, lat2rads, deltaLat, lat2rads, deltaLng,
    a, c, dist_metre, R;

    // Avoid to return NAN, if finding distance between same lat long.
    if(lat1 == lat2 && lon1 == lon2) {
        return 0;
    }

    //Earth Radius (in metre)
    R = 6371000

    lat1rads = this.degreesToRadians(lat1)
    lat2rads = this.degreesToRadians(lat2)
    deltaLat = this.degreesToRadians((lat2-lat1))
    deltaLng = this.degreesToRadians((lng2-lng1))

    a = Math.sin(deltaLat/2) * Math.sin(deltaLat/2) +
        Math.cos(lat1rads) * Math.cos(lat2rads) * Math.sin(deltaLng/2) * Math.sin(deltaLng/2)
    c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a))

    dist_metre = R * c;

    if(isNaN(dist_metre)) {
        return 0;
    }

    return dist_metre
}

LatLong.prototype.degreesToRadians = function (degree) {
    return  degree * Math.PI/180;
}

LatLong.prototype.radiansToDegrees = function(radians) {
    return radians * 180 / Math.PI;
}

/**
 * returns the lat and long of destination point
 * given the start lat, long, aziuth, and distance.
 *
 * @param lat
 * @param lng
 * @param azimuth
 * @param distance
 * @return {*[]}
 */
LatLong.prototype._getDestinationLatLong = function(lat, lng, azimuth, distance_metre)
{
    var lat2, lng2, R, brng, d_km, lat1, lng1 ;

    R = 6378.1 //Radius of the Earth in km

    //Bearing is degrees converted to radians.
    brng = this.degreesToRadians(azimuth);
    d_km = distance_metre/1000;
    lat1 = this.degreesToRadians(lat)
    lng1 = this.degreesToRadians(lng)

    lat2 = Math.asin(Math.sin(lat1) * Math.cos(d_km/R) +
            Math.cos(lat1)* Math.sin(d_km/R)* Math.cos(brng))
    lng2 = lng1 +
            Math.atan2(
                Math.sin(brng) * Math.sin(d_km/R)* Math.cos(lat1),
                Math.cos(d_km/R)- Math.sin(lat1)* Math.sin(lat2));

    //convert back to degrees
    lat2 = this.radiansToDegrees(lat2)
    lng2 = this.radiansToDegrees(lng2)

    return [parseFloat(lat2.toFixed(6)), parseFloat(lng2.toFixed(6))]
}

/**
 * calculates the azimuth in degrees from start point to end point
 *
 * @param lat1
 * @param lng1
 * @param lat2
 * @param lng2
 * @return {*}
 */
LatLong.prototype.calculateBearing = function (lat1,lng1,lat2,lng2)
{
    var startLat, startLong, endLat, endLong, dLong, dPhi, bearing;

    startLat = this.degreesToRadians(lat1)
    startLong = this.degreesToRadians(lng1)
    endLat = this.degreesToRadians(lat2)
    endLong = this.degreesToRadians(lng2)

    dLong = endLong - startLong
    dPhi = Math.log(Math.tan(endLat/2.0+Math.PI/4.0)/Math.tan(startLat/2.0+Math.PI/4.0));

    if (Math.abs(dLong) > Math.PI) {
        if (dLong > 0) {
            dLong = -(2.0 * Math.PI - dLong)
        } else {
            dLong = (2.0 * Math.PI + dLong)
        }
    }

    bearing = (this.radiansToDegrees(Math.atan2(dLong, dPhi)) + 360.0) % 360.0;

    return bearing;
}

/**
 * Invoke to returns every coordinate pair in-between two coordinate pairs
 * given the desired interval
 *
 * @param interval
 * @param azimuth
 * @param lat1
 * @param lng1
 * @param lat2
 * @param lng2
 *
 * @return {Array}
 */
LatLong.prototype._buildCoordinates = function (interval, azimuth, lat1, lng1, lat2, lng2)
{
    var d, dist, counter, coords, range_list, _coord ;

    d = this.getPathLength(lat1, lng1, lat2, lng2);

    dist = parseInt(d/interval);

    coords = [];
    coords.push([lat1,lng1]);
    //:::::::::::::::::::::::::::::::::::::
    //::::: Lodash/UnderScore lib     :::::
    //:::::::::::::::::::::::::::::::::::::
    range_list = _.range(0, dist);

    counter = parseFloat(interval)

    for (var key in range_list) {
        _coord = this._getDestinationLatLong(lat1, lng1, azimuth, counter)
        counter = counter + parseFloat(interval)
        coords.push(_coord);
    }

    coords.push([lat2,lng2])

    return coords;

}

/**
 * Invoke to get coordinates between two location
 * @param lat1
 * @param lng1
 * @param lat2
 * @param lng2
 * @param interval_meters
 * @return {Array|*}
 */
LatLong.prototype.getCoordinates = function (lat1, lng1, lat2, lng2, interval_meters)
{
    var azimuth, coords;

    // point interval in meters
    if (!interval_meters) {
        interval_meters = 20.0
    }

    azimuth = this.calculateBearing(lat1, lng1, lat2, lng2)
    coords = this._buildCoordinates(interval_meters, azimuth, lat1, lng1, lat2, lng2)

    return coords;
}


module.exports = LatLong;
0

In the python code, replace in the function getDestinationLatLong

d = distance/1000

with

d = float(distance)/1000

Or d will be rounded ....

0

Thanks @Peter for the code! I converted it to java. Comparing the results from both, there's a slight difference but they are accurate up to 5 decimal points, which is roughly 0.11 meters. Here's the code:

import java.util.ArrayList;

public class GpsLocations {

    private static final long RADIUS_OF_EARTH = 6371000; // radius of earth in m

    static class MockLocation {
        double lat;
        double lng;

        public MockLocation(double lat, double lng) {
            this.lat = lat;
            this.lng = lng;
        }

        @Override
        public String toString() {
            return "(" + lat + "," + lng + ")";
        }
    }

    public static void main(String args[]) {
        // point interval in meters
        int interval = 1;
        // direction of line in degrees
        //start point
        double lat1 = 43.97076;
        double lng1 = 12.72543;
        // end point
        double lat2 = 43.969730;
        double lng2 = 12.728294;

        MockLocation start = new MockLocation(lat1, lng1);
        MockLocation end = new MockLocation(lat2, lng2);
        double azimuth = calculateBearing(start, end);
        System.out.println(azimuth);
        ArrayList<MockLocation> coords = getLocations(interval, azimuth, start, end);

        for (MockLocation mockLocation : coords) {
            System.out.println(mockLocation.lat + ", " + mockLocation.lng);
        }

    }

    /**
     * returns every coordinate pair in between two coordinate pairs given the desired interval
     * @param interval
     * @param azimuth
     * @param start
     * @param end
     * @return
     */
    private static ArrayList<MockLocation> getLocations(int interval, double azimuth, MockLocation start, MockLocation end) {
        System.out.println("getLocations: " +
                "\ninterval: " + interval +
                "\n azimuth: " + azimuth +
                "\n start: " + start.toString());

        double d = getPathLength(start, end);
        int dist = (int) d / interval;
        int coveredDist = interval;
        ArrayList<MockLocation> coords = new ArrayList<>();
        coords.add(new MockLocation(start.lat, start.lng));
        for(int distance = 0; distance < dist; distance += interval) {
            MockLocation coord = getDestinationLatLng(start.lat, start.lng, azimuth, coveredDist);
            coveredDist += interval;
            coords.add(coord);
        }
        coords.add(new MockLocation(end.lat, end.lng));

        return coords;

    }

    /**
     * calculates the distance between two lat, long coordinate pairs
     * @param start
     * @param end
     * @return
     */
    private static double getPathLength(MockLocation start, MockLocation end) {
        double lat1Rads = Math.toRadians(start.lat);
        double lat2Rads = Math.toRadians(end.lat);
        double deltaLat = Math.toRadians(end.lat - start.lat);

        double deltaLng = Math.toRadians(end.lng - start.lng);
        double a = Math.sin(deltaLat/2) * Math.sin(deltaLat/2) + Math.cos(lat1Rads) * Math.cos(lat2Rads) * Math.sin(deltaLng/2) * Math.sin(deltaLng/2);
        double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
        double d = RADIUS_OF_EARTH * c;
        return d;
    }

    /**
     * returns the lat an long of destination point given the start lat, long, aziuth, and distance
     * @param lat
     * @param lng
     * @param azimuth
     * @param distance
     * @return
     */
    private static MockLocation getDestinationLatLng(double lat, double lng, double azimuth, double distance) {
        double radiusKm = RADIUS_OF_EARTH / 1000; //Radius of the Earth in km
        double brng = Math.toRadians(azimuth); //Bearing is degrees converted to radians.
        double d = distance / 1000; //Distance m converted to km
        double lat1 = Math.toRadians(lat); //Current dd lat point converted to radians
        double lon1 = Math.toRadians(lng); //Current dd long point converted to radians
        double lat2 = Math.asin(Math.sin(lat1) * Math.cos(d / radiusKm) + Math.cos(lat1) * Math.sin(d / radiusKm) * Math.cos(brng));
        double lon2 = lon1 + Math.atan2(Math.sin(brng) * Math.sin(d / radiusKm) * Math.cos(lat1), Math.cos(d / radiusKm) - Math.sin(lat1) * Math.sin(lat2));
        //convert back to degrees
        lat2 = Math.toDegrees(lat2);
        lon2 = Math.toDegrees(lon2);
        return new MockLocation(lat2, lon2);
    }

    /**
     * calculates the azimuth in degrees from start point to end point");
     double startLat = Math.toRadians(start.lat);
     * @param start
     * @param end
     * @return
     */
    private static double calculateBearing(MockLocation start, MockLocation end) {
        double startLat = Math.toRadians(start.lat);
        double startLong = Math.toRadians(start.lng);
        double endLat = Math.toRadians(end.lat);
        double endLong = Math.toRadians(end.lng);
        double dLong = endLong - startLong;
        double dPhi = Math.log(Math.tan((endLat / 2.0) + (Math.PI / 4.0)) / Math.tan((startLat / 2.0) + (Math.PI / 4.0)));
        if (Math.abs(dLong) > Math.PI) {
            if (dLong > 0.0) {
                dLong = -(2.0 * Math.PI - dLong);
            } else {
                dLong = (2.0 * Math.PI + dLong);
            }
        }
        double bearing = (Math.toDegrees(Math.atan2(dLong, dPhi)) + 360.0) % 360.0;
        return bearing;
    }

}
  • Thank you for this code but when I ported this into c# then line drawn using these points are slightly going in the anti-clockwise direction. Please help! – Manish Jain Apr 19 '19 at 11:08
0

Thanks to @PMK here is the C# implementation of it. My change is to get x-Number of points between start and end instead of getting all points in x distance

/// <summary>
        /// Calculates the distance between two lat, long coordinate pairs
        /// </summary>
        /// <param name="lat1"></param>
        /// <param name="lng1"></param>
        /// <param name="lat2"></param>
        /// <param name="lng2"></param>
        /// <returns></returns>
        private double GetPathLength(double lat1, double lng1, double lat2, double lng2)
        {
            var radiusOfEarthInMeters = 6371000;
            var lat1rads = lat1.ToRadian();
            var lat2rads = lat2.ToRadian();
            var deltaLat = (lat2 - lat1).ToRadian();
            var deltaLng = (lng2 - lng1).ToRadian();
            var a = Math.Sin(deltaLat / 2) * Math.Sin(deltaLat / 2) + Math.Cos(lat1rads) *
        Math.Cos(lat2rads) * Math.Sin(deltaLng / 2) * Math.Sin(deltaLng / 2);
            var c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
            var d = radiusOfEarthInMeters * c;

            return d;
        }

        /// <summary>
        /// Returns the lat and long of destination point given the start lat, long, aziuth, and distance
        /// </summary>
        /// <param name="startLat"></param>
        /// <param name="startLong"></param>
        /// <param name="azimuth"></param>
        /// <param name="totalTravelDistance"></param>
        /// <returns></returns>
        private Coordinates GetDestinationLatLong(double startLat,
            double startLong,
            double azimuth,
            double totalTravelDistance)
        {
            var radiusOfEarthInKiloMeters = 6378.1;
            var brng = azimuth.ToRadian(); //Bearing is degrees converted to radians.
            var distanceInKiloMeters = totalTravelDistance / 1000; //Distance m converted to km
            var lat1 = startLat.ToRadian();//Current dd lat point converted to radians
            var lon1 = startLong.ToRadian();//Current dd long point converted to radians
            var lat2 = Math.Asin(Math.Sin(lat1) * Math.Cos(distanceInKiloMeters / radiusOfEarthInKiloMeters)
                + Math.Cos(lat1) * Math.Sin(distanceInKiloMeters / radiusOfEarthInKiloMeters) * Math.Cos(brng));

            var lon2 = lon1 + Math.Atan2(Math.Sin(brng) * Math.Sin(distanceInKiloMeters / radiusOfEarthInKiloMeters) * Math.Cos(lat1),
                   Math.Cos(distanceInKiloMeters / radiusOfEarthInKiloMeters) - Math.Sin(lat1) * Math.Sin(lat2));

            //# convert back to degrees
            lat2 = lat2.ToDegree();
            lon2 = lon2.ToDegree();

            return new Coordinates(lat2, lon2);
        }

        /// <summary>
        /// Calculates the azimuth in degrees from start point to end point
        /// </summary>
        /// <param name="lat1"></param>
        /// <param name="lng1"></param>
        /// <param name="lat2"></param>
        /// <param name="lng2"></param>
        /// <returns></returns>
        public double CalculateBearing(double lat1, double lng1, double lat2, double lng2)
        {
            var startLat = lat1.ToRadian();
            var startLong = lng1.ToRadian();
            var endLat = lat2.ToRadian();
            var endLong = lng2.ToRadian();
            var dLong = endLong - startLong;
            var dPhi = Math.Log(Math.Tan(endLat / 2.0 + Math.PI / 4.0) / Math.Tan(startLat / 2.0 + Math.PI / 4.0));
            if (Math.Abs(dLong) > Math.PI)
            {
                if (dLong > 0.0)
                {
                    dLong = -(2.0 * Math.PI - dLong);
                }
                else
                {
                    dLong = (2.0 * Math.PI + dLong);
                }
            }
            var bearing = (Math.Atan2(dLong, dPhi).ToDegree() + 360.0) % 360.0;
            return bearing;
        }

Here is how I call it,

//params lat1, lng1, lat2, lng2

 var totalRouteDistance = GetPathLength(lat1, lng1, lat2, lng2);
            //var dist = (totalRouteDistance % interval);
            var dist = (totalRouteDistance / noOfPoints);
            var interval = Math.Round(dist, 0);
            var counter = interval;
            var coords = new List<Coordinates>();
           // coords List object
            foreach (var distance in Enumerable.Range(0, noOfPoints))
            {
                var coord = GetDestinationLatLong(lat1, lng1, azimuth, counter);
                counter = counter + interval;
                coords.Add(coord);
            }

Not the answer you're looking for? Browse other questions tagged or ask your own question.