# What Python package to use to calculate midpoint?

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

• arcpy from ArcGIS software installation, but I'm assuming your are looking for something open-source. There's probably a module from QGIS if that helps. Mar 22, 2013 at 17:33

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

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.

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."""
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))
``````
• Note that this library depends and uses GeographicLib internally, so they should get the same results Mar 15, 2016 at 2:22
• @MikeT True about the dependency, but from my understanding it does the math differently, using 3-dimensional surface normal vectors instead of spherical lat/lon coordinates, and that's how it avoids problems at the poles. they should get the same results, too :) Mar 15, 2016 at 5:53

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.

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)