使用python改进多边形计算中的点

Question

对于我在 python 中的项目，我需要定义大量坐标（纬度、经度）属于哪个国家。我已经设法使用 shapely Point.within(Polygon) 方法（请参见下面的代码）这样​​做，提供了从 natural earth 下载的国家边界的 shapefile（必须首先将 shapefile 拆分为单部分多边形，我没有找到如何正确处理多部分多边形）。

虽然该方法运行良好，但进行大量查询时有点慢。性能可能与 shapefile 分辨率密切相关，但需要精度。我已经在边界框预选回合中取得了一些进展（只有在边界框内具有坐标的多边形被选中），但我正在尝试进一步改进它。我该如何进行？我正在考虑使用内部边界框来快速区分明确的点，但后来我不知道如何构建它们。或者也许最好一劳永逸地制作一些花哨的查找表？我必须说我不熟悉哈希表和类似的东西，这是一个选择吗？

谢谢

#env python3
from shapely.geometry import Polygon, Point
import shapefile

class CountryLoc:
    def __init__(self, shppath, alpha3pos = 1):
        polygon = shapefile.Reader(shppath)
        self.sbbox = [pol.bbox for pol in polygon.shapes()] 
        self.shapePoints = [pol.points for pol in polygon.shapes()]
        self.shapeCnames = [pol.record[alpha3pos] for pol in polygon.shapeRecords()] 
        self.shapecount = polygon.numRecords

    def which_country(self,lat,lon): 
        countries = self._pointcountry(lat,lon)
        if len(countries) > 1:
            print('too many countries at coord ({:.3f},{:.3f}) : {}'.format(lat,lon,countries)) 
        try:
            return countries[0]
        except:
            #print('no country ref @ latlon ({},{})'.format(lat,lon))
            return 'OXO'

    def _findBboxMatch(self,lat,lon):
        matchidslist = []
        for ids in range(self.shapecount):
            if lat >= self.sbbox[ids][1] and lat <= self.sbbox[ids][3]:
                if self.sbbox[ids][0] > self.sbbox[ids][2] :
                # rem RUS and USA BBOX are spanning accross the +180-180 cut 
                    if not(lon >= self.sbbox[ids][2] and lon <= self.sbbox[ids][0]):    
                        matchidslist.append(ids)            
                else:       
                    if lon >= self.sbbox[ids][0] and lon <= self.sbbox[ids][2]: 
                        matchidslist.append(ids)
        return matchidslist

    def _pointcountry(self,lat,lon):
        coord = Point(lon,lat)
        matchidlist = self._findBboxMatch(lat,lon) ## BBOX PRESELECT
        matchcountrylist = []
        for ids in matchidlist:
            pp = Polygon(self.shapePoints[ids])
            if coord.within(pp):
                matchcountrylist.append(self.shapeCnames[ids])
        return matchcountrylist

    def printCountry(self,lat,lon): 
        ctry = self.which_country(lat,lon)
        print('coord. ({},{}) is in {}'.format(lat,lon,ctry))
        bbmatch = self._findBboxMatch(lat,lon)
        print('matching BBOXs are {}'.format([self.shapeCnames[k] for k in bbmatch]))
        print(' ')

if __name__ == '__main__':
    # testing
    cc = input('lat,lon ? ')
    coords = [float(cc.split(',')[0]) , float(cc.split(',')[1])]
    coloc = CountryLoc('./borders_shapefile.shp', alpha3pos=9)
    coloc.printCountry(coords[0],coords[1])

Answer 1

您需要一个加速结构，例如quadtree、k-d tree、spatial hash table 等。

设置形状数据时，您将根据形状在平面中的位置加载结构。例如，对于四叉树，您递归地将空间细分为四个象限，直到每个叶象限与少量形状重叠（或它们都不重叠）。然后在每个叶子上记录一个形状引用列表。

稍后，当您搜索与特定点重叠的形状时，您会根据每个大约 log n 级别的两个比较操作遍历细分树。当您到达正确的叶子象限时，您只需使用Point.within 函数检查少量形状。