9

我正在寻找 python 的地理库。我需要能够执行以下操作：

  1. 使用大圆距离（不是线性距离计算）获取 2 点之间的距离（以米为单位）
  2. 检查点是否在多边形内
  3. 每秒执行 1 和 2 数千次

一开始我看过这篇文章： Python module for storage and querying geocoses and start to use geopy。我遇到了2个问题：

  1. Geopy 不支持多边形
  2. geoPy的高CPU使用率（计算一个点和相对5000个点之间的距离大约需要140ms的CPU）

我一直在寻找并找到最好的 Python GIS 库？和https://gis.stackexchange.com/。它看起来很有希望，因为 geos 正在使用经过编译的 C 代码，它应该更快并且支持多边形。问题是 geos/OGR 执行线性距离计算而不是球体。这消除了所有其他基于地理的模块（如 GEODjango 和 shapely）。我在这里错过了什么吗？我不认为我是第一个使用 python 执行 GIS 计算并希望得到准确结果的人。

4

2 回答 2

1

更新

现在继续完成该库中的其他 576 个函数，不包括已完成的两个多边形函数、已完成的三个球体距离算法，以及两个新的，一个 angle_box_2d 和 angle_contains_ray_2d。另外，我切换到 C 版本，这样就不需要 externs，简化了工作。将旧的 C++ 版本放在 old_c++ 目录中，所以它仍然存在。

测试的性能，它与答案底部列出的相同。

更新 2

所以只是一个快速更新，我还没有完成整个库（我只完成了大约 15%），但我已经添加了这些未经测试的函数，以防你立即需要它们，在 github 上，添加到多边形和球体距离算法中的旧点。

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

我在上面评论过的那些可能不起作用，其他可能不起作用，但同样 - 多边形和球体距离肯定会起作用。您可以指定米、公里、英里、海里，这对球面距离无关紧要，输出与输入的单位相同——算法与单位无关。

我今天早上把它放在一起，所以它目前只提供多边形中的点、凸多边形中的点和三种不同类型的球面距离算法，但至少你要求的那些现在在那里供你使用。我不知道是否与那里的任何其他 python 库存在名称冲突，这些天我只在外围参与 python，所以如果有更好的名称，我愿意接受建议。

在github上：https ://github.com/hoonto/pygeometry

它只是此处描述和实现的功能的 python 桥梁：

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

GEOMETRY 库实际上非常好，所以我认为为 python 连接所有这些函数会很有用，我今晚可能会这样做。

编辑：其他几件事

  1. 因为数学函数实际上是编译的 C++，所以你当然需要确保共享库在路径中。您可以修改geometry.py 以指向您想要放置该共享库的任何位置。
  2. 仅针对 linux 编译，.o 和 .so 在 x86_64 fedora 上编译。
  3. 球面距离算法需要弧度，因此您需要将十进制纬度/经度度数转换为弧度，如geometry.py 中所示。

如果您在 Windows 上确实需要此功能，请告诉我，只需几分钟即可在 Visual Studio 中完成。但除非有人问，否则我可能暂时不理会它。

希望这可以帮助！

Rgds....Hoonto/马特

（新提交：SHA：4fa2dbbe849c09252c7bd931edfe8db478de28e6 - 修复了一些问题，如弧度转换以及 py 函数的返回类型。还添加了一些基本性能测试以确保库正常运行。）

测试结果 在每次迭代中，一次调用 sphere_distance1 一次调用 polygon_contains_point_2d，因此总共调用了 2 次库。

  • ~0.062s：2000 次迭代，4000 次调用
  • ~0.603s：20000 次迭代，40000 次调用
  • ~0.905s：30000 次迭代，60000 次调用
  • ~1.198s：40000 次迭代，80000 次调用
于 2013-06-18T17:53:15.580 回答
0

如果球面计算就足够了，我只会使用 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)。

于 2013-06-19T10:57:06.157 回答