使用 python 改进多边形计算中的点
Improving point within polygon calculations with python
对于我在 python 的项目,我需要定义大量坐标 (lat,lon) 属于哪个国家/地区。我已经设法使用 shapely Point.within(Polygon) 方法(见下面的代码)做到这一点,提供了从 natural earth 下载的国家边界的 shapefile(shapefile 必须首先拆分为单个部分的多边形,我没有'找到如何正确处理多部分多边形)。
虽然该方法工作正常,但进行大量查询时有点慢。性能可能与 shapefile 分辨率密切相关,但需要精度。我已经在边界框预选回合中取得了一些进展(仅检查其边界框内具有坐标的多边形),但我正在尝试进一步改进它。我该如何进行?我正在考虑使用内部边界框来快速区分明确的点,但后来我不知道如何构建它们。或者也许最好一劳永逸地进行一些花哨的查找 table ?我必须说我不熟悉 hash tables 之类的,它是一个选项吗?
谢谢
#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])
你想要一个加速结构,比如quadtree, a k-d tree, a spatial hash table,等等。
设置形状数据时,您将根据形状在平面中的位置加载结构。例如,对于四叉树,您递归地将 space 细分为四个象限,直到每个叶象限与少量形状(或其中的 none 重叠)。然后你在每个叶子上记录一个形状参考列表。
稍后,当您搜索与特定点重叠的形状时,您可以根据每个大约 log n 级别的两个比较操作来遍历细分树.当你到达正确的叶子象限时,你只需要用Point.within
函数检查少量的形状。
对于我在 python 的项目,我需要定义大量坐标 (lat,lon) 属于哪个国家/地区。我已经设法使用 shapely Point.within(Polygon) 方法(见下面的代码)做到这一点,提供了从 natural earth 下载的国家边界的 shapefile(shapefile 必须首先拆分为单个部分的多边形,我没有'找到如何正确处理多部分多边形)。
虽然该方法工作正常,但进行大量查询时有点慢。性能可能与 shapefile 分辨率密切相关,但需要精度。我已经在边界框预选回合中取得了一些进展(仅检查其边界框内具有坐标的多边形),但我正在尝试进一步改进它。我该如何进行?我正在考虑使用内部边界框来快速区分明确的点,但后来我不知道如何构建它们。或者也许最好一劳永逸地进行一些花哨的查找 table ?我必须说我不熟悉 hash tables 之类的,它是一个选项吗?
谢谢
#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])
你想要一个加速结构,比如quadtree, a k-d tree, a spatial hash table,等等。
设置形状数据时,您将根据形状在平面中的位置加载结构。例如,对于四叉树,您递归地将 space 细分为四个象限,直到每个叶象限与少量形状(或其中的 none 重叠)。然后你在每个叶子上记录一个形状参考列表。
稍后,当您搜索与特定点重叠的形状时,您可以根据每个大约 log n 级别的两个比较操作来遍历细分树.当你到达正确的叶子象限时,你只需要用Point.within
函数检查少量的形状。