Fast Python Gis Library That Supports Great Circle Distance And Polygon

I was looking for a geographical library for python. I need to be able to do the following: Get the distance between 2 points (in meters) using Great-circle distance (not liner di

Solution 1:

UPDATE

Moving on now to finishing out the other 576 functions in that library not including the two polygon functions that are finished, the three sphere distance algorithms that are done, and two new ones, an angle_box_2d and angle_contains_ray_2d. Also, I switched to the C version so that externs are not needed, simplifies the work. Put the old C++ version in directory old_c++, so its still there.

Tested performance, it is identical as listed at the bottom of the answer.

---

UPDATE 2

So just a quick update, I haven't finished the whole library yet (I'm only about 15% of the way through), but I've added these untested functions, in case you need them right away, on github, to add to the old point in polygon and sphere distance algorithms.

angle_box_2d
angle_contains_ray_2d
angle_deg_2d
angle_half_2d # MLM: double *
angle_rad_2d
angle_rad_3d
angle_rad_nd
angle_turn_2d
anglei_deg_2d
anglei_rad_2d
annulus_area_2d
annulus_sector_area_2d
annulus_sector_centroid_2d # MLM: double *
ball_unit_sample_2d # MLM: double *
ball_unit_sample_3d # MLM: double *
ball_unit_sample_nd # MLM; double *
basis_map_3d #double *
box_01_contains_point_2d
box_01_contains_point_nd
box_contains_point_2d
box_contains_point_nd
box_ray_int_2d
box_segment_clip_2d
circle_arc_point_near_2d
circle_area_2d
circle_dia2imp_2d
circle_exp_contains_point_2d
circle_exp2imp_2d
circle_imp_contains_point_2d
circle_imp_line_par_int_2d
circle_imp_point_dist_2d
circle_imp_point_dist_signed_2d
circle_imp_point_near_2d
circle_imp_points_2d # MlM: double *
circle_imp_points_3d # MLM: double *
circle_imp_points_arc_2d
circle_imp_print_2d
circle_imp_print_3d
circle_imp2exp_2d
circle_llr2imp_2d # MLM: double *
circle_lune_area_2d
circle_lune_centroid_2d # MLM; double *
circle_pppr2imp_3d

The ones that I've commented above probably won't work, the others might, but again - polygon & sphere distances definitely do. And you can specify meters, kilometers, miles, nautical miles, it doesn't really matter on the spherical distance ones, the output is in the same units as the input - the algorithms are agnnostic to the units.

---

I put this together this morning so it currently only provides the point in polygon, point in convex polygon, and three different types of spherical distance algorithms, but at least those ones that you requested are there for you to use now. I don't know if there is a name conflict with any other python library out there, I only get peripherally involved with python these days, so if there's a better name for it I'm open to suggestions.

On github: https://github.com/hoonto/pygeometry

It is just a python bridge to the functions described and implemented here:

http://people.sc.fsu.edu/~jburkardt/cpp_src/geometry/geometry.html

The GEOMETRY library is pretty good actually, so I think it'll be useful to bridge all of those functions for python, which I'll do probably tonight.

Edit: a couple other things

1. Because the math functions are actually compiled C++, you do of course need to make sure that the shared library is in the path. You can modify the geometry.py to point at wherever you want to put that shared library though.
2. Only compiled for linux, the .o and .so were compiled on x86_64 fedora.
3. The spherical distance algorithms expect radians so you need to convert decimal lat/lon degrees for example to radians, as shown in geometry.py.

If you do need this on Windows let me know, it should only take a couple minutes to get it worked out in Visual Studio. But unless someone asks I'll probably just leave it alone for now.

Hope this helps!

Rgds....Hoonto/Matt

(new commit: SHA: 4fa2dbbe849c09252c7bd931edfe8db478de28e6 - fixed some things, like radian conversions and also the return types for the py functions. Also added some basic performance tests to make sure the library performs appropriately.)

Test Results In each iteration, one call to sphere_distance1 and one call polygon_contains_point_2d so 2 calls to the library total.

• ~0.062s : 2000 iterations, 4000 calls
• ~0.603s : 20000 iterations, 40000 calls
• ~0.905s : 30000 iterations, 60000 calls
• ~1.198s : 40000 iterations, 80000 calls

Solution 2:

If spherical calculation is enough I'd just use numpy for distance and matplotlib for polygon check (as you find similar proposals in stackoverflow).

from math import asin, cos, radians, sin, sqrt
import numpy as np

defgreat_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) ** 2return2 * asin(min(1, sqrt(a))) * radius

defgreat_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

defaux_generateLatLon():
    import random
    while1:
        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 inrange(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 performanceimport 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 distanceassert(abs(dMin1 - dMin2) < 0.0001)
    printprint"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 Polesfrom 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))

If you want to do more, the Quantum GIS toolset probably is worth a look: PyQGIS Developer Cookbook (docs.qgis.org).