球状の計算で十分な場合は、距離には numpy を使用し、多角形チェックには matplotlib を使用します (stackoverflow で同様の提案があるため)。
from math import asin, cos, radians, sin, sqrt
import numpy as np
def great_circle_distance_py(pnt1, pnt2, radius):
""" Returns distance on sphere between points given as (latitude, longitude) in degrees. """
lat1 = radians(pnt1[0])
lat2 = radians(pnt2[0])
dLat = lat2 - lat1
dLon = radians(pnt2[1]) - radians(pnt1[1])
a = sin(dLat / 2.0) ** 2 + cos(lat1) * cos(lat2) * sin(dLon / 2.0) ** 2
return 2 * asin(min(1, sqrt(a))) * radius
def great_circle_distance_numpy(pnt1, l_pnt2, radius):
""" Similar to great_circle_distance_py(), but working on list of pnt2 and returning minimum. """
dLat = np.radians(l_pnt2[:, 0]) - radians(pnt1[0]) # slice latitude from list of (lat, lon) points
dLon = np.radians(l_pnt2[:, 1]) - radians(pnt1[1])
a = np.square(np.sin(dLat / 2.0)) + np.cos(radians(pnt1[0])) * np.cos(np.radians(l_pnt2[:, 0])) * np.square(np.sin(dLon / 2.0))
return np.min(2 * np.arcsin(np.minimum(np.sqrt(a), len(a)))) * radius
def aux_generateLatLon():
import random
while 1:
yield (90.0 - 180.0 * random.random(), 180.0 - 360.0 * random.random())
if __name__ == "__main__":
## 1. Great-circle distance
earth_radius_m = 6371000.785 # sphere of same volume
nPoints = 1000
nRep = 100 # just to measure time
# generate a point and a list of to check against
pnt1 = next(aux_generateLatLon())
l_pnt2 = np.array([next(aux_generateLatLon()) for i in range(nPoints)])
dMin1 = min([great_circle_distance_py(pnt1, pnt2, earth_radius_m) for pnt2 in l_pnt2])
dMin2 = great_circle_distance_numpy(pnt1, l_pnt2, earth_radius_m)
# check performance
import timeit
print "random points: %7i" % nPoints
print "repetitions : %7i" % nRep
print "function 1 : %14.6f s" % (timeit.timeit('min([great_circle_distance_py(pnt1, pnt2, earth_radius_m) for pnt2 in l_pnt2])', 'from __main__ import great_circle_distance_py , pnt1, l_pnt2, earth_radius_m', number=nRep))
print "function 2 : %14.6f s" % (timeit.timeit('great_circle_distance_numpy(pnt1, l_pnt2, earth_radius_m)' , 'from __main__ import great_circle_distance_numpy, pnt1, l_pnt2, earth_radius_m', number=nRep))
# tell distance
assert(abs(dMin1 - dMin2) < 0.0001)
print
print "min. distance: %14.6f m" % dMin1
## 2. Inside polygon?
# Note, not handled:
# - the "pathological case" mentioned on http://paulbourke.net/geometry/polygonmesh/
# - special situations on a sphere: polygons covering "180 degrees longitude edge" or the Poles
from matplotlib.path import Path
x = y = 1.0
l_pnt2 = [(-x, -y), (x, -y), (x, y), (-x, y), (-x, -y)]
path = Path(l_pnt2)
print "isInside ?"
for pnt in [(0.9, -1.9), (0.9, -0.9)]:
print " ", pnt, bool(path.contains_point(pnt))
もっとやりたい場合は、Quantum GIS ツールセットは一見の価値があります: PyQGIS Developer Cookbook (docs.qgis.org)。