# Getting all vertex lat, long coordinates every 1 meter between two known points using PHP?

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
• 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

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
a = sin(deltaLat/2) * sin(deltaLat/2) + cos(lat1rads) *
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
d = distance/1000 #Distance m converted to km

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
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
d = distance/1000 #Distance m converted to km
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'''
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

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) && (!isset(\$pregResultB)
|| strstr(\$preResultA,\$pregResultB) == 0
|| !\$pointset))) {
\$numberString = \$pregResultA;
\$numberString = str_replace('.','',\$numberString);
\$numberString = str_replace(',','.',\$numberString);
}
elseif (isset(\$pregResultB) && (!isset(\$pregResultA)
|| strstr(\$pregResultB,\$preResultA) == 0
|| !\$commaset)) {
\$numberString = \$pregResultB;
\$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
\$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
\$d = \$distance / 1000; #Distance m converted to km
\$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

return [\$lat2, \$lon2];
}

function calculateBearing(\$lat1,\$lng1,\$lat2,\$lng2){
// '''calculates the azimuth in degrees from start point to end point'''
\$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;
\$dist = \$rapydUnpack;
\$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

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)
{
a, c, dist_metre, R;

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

R = 6371000

a = Math.sin(deltaLat/2) * Math.sin(deltaLat/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
}

return  degree * Math.PI/180;
}

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.
d_km = distance_metre/1000;

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

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;

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;
``````

In the python code, replace in the function getDestinationLatLong

``````d = distance/1000
``````

with

``````d = float(distance)/1000
``````

Or d will be rounded ....

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<>();
for(int distance = 0; distance < dist; distance += interval) {
MockLocation coord = getDestinationLatLng(start.lat, start.lng, azimuth, coveredDist);
coveredDist += interval;
}

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 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 d = distance / 1000; //Distance m converted to km
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");
* @param start
* @param end
* @return
*/
private static double calculateBearing(MockLocation start, MockLocation end) {
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

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 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 distanceInKiloMeters = totalTravelDistance / 1000; //Distance m converted to km
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 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;