Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

share|improve this question
    
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

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.

share|improve this answer

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))
share|improve this answer
1  
Note that this library depends and uses GeographicLib internally, so they should get the same results – Mike T Mar 15 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 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.

share|improve this answer

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()
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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