[GIS] How to Specify Haversine when using Buffer Method in Shapely and how to get Haversine distance between two Shapely Point objects

Question

I have a series of latitude/longitude pairs and I form Point() objects for each pair. I would like to form buffer regions around each point using a Haversine distance argument instead of Euclidean distance. That way, points on the boundary of the buffer are an equal Haversine distance away from the center of the buffer. Is there a way I can do this?

Also, is there a built-in method for getting the Haversine distance between two Shapely Point() objects? I know I can use the Haversine formula for the lat/long pairs, but I was wondering if there was anything built-in.

Answer 1

Instead of Haversine distance formula, you can use directly pyproj python module. For instance, for producing a 1000 m buffer (EPSG:32612) around this point (-112.171129818, 40.1731054115), EPSG:4326, next code does the job.

from shapely.geometry import Point, Polygon

import pyproj

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

pt = (-112.171129818, 40.1731054115)

x, y = pyproj.transform(srcProj, dstProj, pt[0], pt[1])

pt = Point(x,y)

buffer = pt.buffer(1000)

buffer_points =  zip(buffer.exterior.coords.xy[0], buffer.exterior.coords.xy[1])

proj_buffer_points = []

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

print Polygon(proj_buffer_points)

After running the code, at console is printed buffer (WKT format) in long/lat coordinates (EPSG:4326). With the help of QuickWKT plugin, it is visualized (EPSG:4326) at next image with QGIS.

When 'on the fly' CRS transformation to EPSG:32612 is enabled, it can be observed circular buffer shape.