I asked a question earlier about how to calculate a midpoint. I was wondering if there is a Python library that does this alreay?
This python package
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
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'])
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.""" # 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))
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) point_object.AddPoint(xMid,yMid) featObj.SetGeometry(point_object) myLayer.CreateFeature(featObj) #close both shapefiles mySource.Destroy() dataSource.Destroy()