[GIS] Getting all vertex lat, long coordinates every 1 meter between two known points using PHP

coordinateslatitude longitudePHP

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.

Best Answer

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.

Related Solutions

[GIS] Converting local Coordinates to lat/long Coordinates using QGIS

You have to build the proj.4 string for the CRS from the parameters given:

+proj=lcc +lat_1=44.883333 +lat_2=45.133333 +lat_0=44.791111 +lon_0=-93.383333 +x_0=152400.000000 +y_0=30480.000000 +a=6378418.9409999996 +b=6357033.3098455509 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

Note that x_0 and y_0 have to be in meters, while false Easting and Northing from the definition (and the WKT mkennedy added) are in units of the projection (us-ft).

With that, you can load the data as delimited text, with X set to the itm column and Y to the utm column. If you have not set QGIS to prompt for the CRS, rightclick on the layer -> Set CRS for layer

to get the right alignment:

enter image description here

I'm not sure why they have named the second column utm, the projection has nothing to do with Universal Tramnsverse Mercator.

Python Haversine – Compute Distance Between Coordinates and Print CSV

If I've understood correctly, this code should do what you're after. Using your input

Distance,Longitude,Latitude
,-77.60483,40.31117
,-77.58167,40.3045
,-77.52883,40.24917
,-77.508,40.14917
,-77.49617,40.13117

I'm producing this output. I've assumed that you want the distances to go forward - i.e. that line 2 has its distance set to the distance to the next point (line 3).The last line is set to 0 as there is no following line to measure to

Distance,Longitude,Latitude
2.0992,-77.6048,40.3112
7.6122,-77.5817,40.3045
11.2593,-77.5288,40.2492
2.2399,-77.5080,40.1492
0.0000,-77.4962,40.1312

Keep your existing code, and append this snippet. The csv library makes it nice and easy to work with CSV files.

I'm also using a deque. This is a FIFO (first in, first out) data structure. I use this to keep a queue of points. Once this has two entries, I compute the distance between the two points, write it out, then discard the oldest point.

import csv
import collections
filein = open("test.csv","rU") # universal newline mode
fileout = open("output.csv","w")
fileout.write("Distance,Longitude,Latitude\n")
csvrdr = csv.reader(filein, delimiter=",")
linecount = 0
# use a deque. once this hits 2 entries, compute the distance
# and then remove the oldest entry
todo = collections.deque()
for line in csvrdr:
    if linecount > 0: # skip header
        _, longitude, latitude = line # _ means discard
        longitude = float(longitude)
        latitude = float(latitude)
        todo.append((longitude, latitude))
    if len(todo)==2:
        # got two points, compute distance and remove the oldest
        oldlon, oldlat = todo.popleft()
        distance = distance_on_unit_sphere(latitude, longitude, oldlat, oldlon)
        fileout.write("%2.4f,%2.4f,%2.4f\n" % (distance, oldlon, oldlat))
    linecount += 1
# got one entry in deque, spit that out with 0
assert(len(todo)==1)
oldlon, oldlat = todo.popleft()
fileout.write("%2.4f,%2.4f,%2.4f\n" % (0.0, oldlon, oldlat))
fileout.close()
filein.close()

Best Answer

Related Solutions

[GIS] Converting local Coordinates to lat/long Coordinates using QGIS

Python Haversine – Compute Distance Between Coordinates and Print CSV

Related Question