What Python package to use to calculate midpoint?

Question

I asked a question earlier about how to calculate a midpoint. I was wondering if there is a Python library that does this alreay?

Answer 1

This python package

http://pypi.python.org/pypi/geographiclib

solves the direct and inverse geodesic problems for an ellipsoid. In your case, solve the inverse problem to get a distance and an azimuth at the first point and then solve the direct problem from the first point with the computed azimuth and half the distance. For a quick illustration of using this package, see

http://geographiclib.sourceforge.net/html/other.html#python

Here, for example, is the calculation of the point midway between JFK and SIN airports

# Compute point midway between JFK and SIN airports
import sys
sys.path.append("/usr/local/lib/python/site-packages")
from geographiclib.geodesic import Geodesic

# Coordinates of airports
lat1,lon1 = 40.640,-73.779 # JFK
lat2,lon2 =  1.359,103.989 # SIN

# Compute path from 1 to 2
g = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2);

# Compute midpoint starting at 1
h1 = Geodesic.WGS84.Direct(lat1, lon1, g['azi1'], g['s12']/2);
print(h1['lat2'],h1['lon2']);

# Alternatively, compute midpoint starting at 2
h2 = Geodesic.WGS84.Direct(lat2, lon2, g['azi2'], -g['s12']/2);
print(h2['lat2'],h2['lon2']);

On my system, this gives

(70.34117458722292, 97.02347775257729)
(70.34117458722295, 97.02347775257725)

ADDENDUM

This calculation is a little simpler with Version 1.46 of the GeographicLib python package

# Define the path from 1 to 2
l = Geodesic.WGS84.InverseLine(lat1, lon1, lat2, lon2)

# Compute the midpoint
m = l.Position(0.5 * l.s13)
print(m['lat2'],m['lon2'])

This gives

(70.34117458722294, 97.02347775257721)

For documentation see http://geographiclib.sourceforge.net/html/python.

Answer 2

I've just come across nvector (http://www.navlab.net/nvector), which seems to do this nicely, including taking care of issues at the poles and the 180/-180 meridian. My code using their methods:

import nvector as nv

class Vertex:
    def __init__( self, lat, lon ):
        self.lat = lat
        self.lon = lon

def GetGeodeticMidpoint(vert1, vert2):
    """Given two Vertices, return the geodetic midpoint of the great circle arc between them, on the WGS84 ellipsoid. Uses nvector."""
    # see http://nvector.readthedocs.org/en/latest/src/overview.html?highlight=midpoint#description    
    wgs84 = nv.FrameE(name='WGS84')
    n_EB_E_t0 = wgs84.GeoPoint(vert1.lat, vert1.lon, degrees=True).to_nvector()
    n_EB_E_t1 = wgs84.GeoPoint(vert2.lat, vert2.lon, degrees=True).to_nvector()
    path = nv.GeoPath(n_EB_E_t0, n_EB_E_t1)    
    halfway = 0.5
    g_EB_E_ti = path.interpolate(halfway).to_geo_point()
    lat_ti, lon_ti = g_EB_E_ti.latitude_deg, g_EB_E_ti.longitude_deg
    return Vertex(float(lat_ti), float(lon_ti))

Answer 3

You could try the Angles library, which does bearing and distance between points on the unit sphere, but I think you'd still have to calculate the midpoint yourself.

Answer 4

GDAL/OGR (Geospatial Data Abstraction Library) has a centroid function.

The following script reads in a line shapefile and writes the midpoint out to a point shapefile (source):

#script that finds the midpoint of transect.shp
#script then creates a point shapefile from the line midpoint

from osgeo import ogr
myDriver = ogr.GetDriverByName('ESRI Shapefile')

#create layer from transect shapefile
dataSource = myDriver.Open('transect.shp',1)
lineLayer = dataSource.GetLayer()

#get geometry from first line feature in layer:
lineFeature = lineLayer.GetFeature(0)
lineGeom = lineFeature.geometry()

#get midpoint of line and get x/y coords
midPoint = lineGeom.Centroid()
xMid = midPoint.GetX()
yMid = midPoint.GetY()

#Create a point shapefile for the midpoint
mySource = myDriver.CreateDataSource('MidPoint.shp')
myLayer = mySource.CreateLayer('layer1',geom_type=ogr.wkbPoint)
featObj = ogr.Feature(myLayer.GetLayerDefn())
point_object = ogr.Geometry(ogr.wkbPoint)
point_object.AddPoint(xMid,yMid)
featObj.SetGeometry(point_object)
myLayer.CreateFeature(featObj)

#close both shapefiles
mySource.Destroy()
dataSource.Destroy()