[GIS] Working with geographic coordinates in shapely

coordinate systempythonshapely

I have a very basic question about working with the shapely module in Python. I have some geographic coordinates received from a GPS. I would like to use shapely to calculate the great cicle distance in meters between two points. I started with:

>>> from shapely.geometry import Point
>>> p1 = Point(43.374880, -78.119956)
>>> p2 = Point(43.374868, -78.119666)

I believe that this gives me two points in the cartesian coordinate system, which isn't going to be very useful. These coordinates come from a GPS which presumably is using the WGS84 CRS. I have seen some example using pyproj and shapely.ops.transform to transform between coordinate systems, but these all seem to involve data points that are already in some valid cartographic coordinate system with an EPSG identifier. I'm not sure what do with naive points like the above that don't have an associated geographic CRS.

How do associate these points with the WGS84 CRS such that p1.distance(p2) will yield the distance in meters? Am I using the wrong tools?

Best Answer

The problem you're trying to solve is the Inverse Geodesic Problem. Happily the Python geographiclib can do this for you:

from geographiclib.geodesic import Geodesic

p1_lat, p1_lon = 43.374880, -78.119956
p2_lat, p2_lon = 43.374868, -78.119666
geod = Geodesic.WGS84

# note the Inverse method expects latitude then longitude
g = geod.Inverse(p1_lat, p1_lon, p2_lat, p2_lon)

print("Distance is {:.2f}m".format(g['s12']))

Giving the result:

Distance is 23.54m

That said, given your points are likely to be only 10s of metres apart, then working in projected units will likely be much easier.

Related Solutions

[GIS] Exporting geographic coordinates in QGis as projected (flat) coordinates

Kml has always degrees as units, so you will have no luck with that format.

Exporting as shapefile should work. Note that you have to choose the CRS again (QGIS does not remember it from the last export).

Can you add the shapefile to the canvas, and look up the extent under properties, metadata tab? Furthermore, you can open the .prj file that was created along with the shapefile. It should contain the projection data for NAD83(HARN) / South Carolina (ft).

[GIS] Shapely python spatial join – point with line

If you have datasets with (long,lat) coordinates, you need to project them prior to distance calculations. For instance, next code opens a point and a line layers previous to these calculations. Results (all distances and minimal distances and its points) are printed by console.

import fiona
from shapely.geometry import shape, Point, LineString
import pyproj

srcProj = pyproj.Proj(init='EPSG:4326')
dstProj = pyproj.Proj(init='EPSG:32612')

path1 = '/home/zeito/pyqgis_data/points_test_4326.shp'
path2 = '/home/zeito/pyqgis_data/new_lines_4326.shp'

points = fiona.open(path1)
line = fiona.open(path2)

points = [ (shape(feat["geometry"]).xy[0][0], shape(feat["geometry"]).xy[1][0]) 
           for feat in points ]

lines = [ zip(shape(feat["geometry"]).coords.xy[0], shape(feat["geometry"]).coords.xy[1]) 
          for feat in line ]

proj_lines = [ [] for i in range(len(lines)) ]

for i, item in enumerate(lines):
    for element in item:
        x = element[0]
        y = element[1]
        x, y = pyproj.transform(srcProj, dstProj, x, y)
        proj_lines[i].append((x, y))

proj_points = []

for point in points:
    x = point[0]
    y = point[1]
    x, y = pyproj.transform(srcProj, dstProj, x, y)    
    proj_points.append(Point(x,y))

distances = [ [] for i in range(len(lines)) ]

for i, line in enumerate(proj_lines):
    for point in proj_points:
        distances[i].append(LineString(line).distance(point))

print distances

for i, list in enumerate(distances):
    print "minimal distance: {:.2f} to this point: {}".format(min(list), points[i])

After running the code, it was gotten:

[[240.76129848041424, 82.17367496992877, 534.8119728316819], [180.2963692590486, 601.7275241365959, 1155.1004277619784]]
minimal distance: 82.17 to this point: (-112.17112916813936, 40.17311785635686)
minimal distance: 180.30 to this point: (-112.1560813471297, 40.170581085323384)

Layers look like at next image. QuickWKT plugin of QGIS helped to corroborate that these determined points are corresponding to minimal distance.

Related Question