4

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. – GIS Danny Mar 22 '13 at 17:33
10

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.

4

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))
  • 1
    Note that this library depends and uses GeographicLib internally, so they should get the same results – Mike T Mar 15 '16 at 2:22
  • 1
    @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 :) – Paulo Raposo Mar 15 '16 at 5:53
0

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.

0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.