現在、GeoJson ファイルがあり、次の関数で shapely を使用しています。
座標を取り、近隣の名前を返します
def get_neighb(lat, lon):
"""Input Latitude and Longitude, Returns Neighborhood Name"""
point = Point(lon, lat)
found = False
for feature in geo_data['features']:
polygon = shape(feature['geometry'])
if polygon.contains(point):
return(feature['properties']['neighborhood'])
found = True
if found is False:
return('NA')
# Initialize list
tn = ['']*data.shape[0]
for i in range(len(tn)):
tn[i] = get_neighb(data.latitude[i], data.longitude[i])
これは機能しますが、非常に遅いため、どのように高速化できるかを考え、現在 4,000,000 行で実行しています。
たとえば、PostGIS データベースの高負荷な機械を避けたい場合は、rtreeパッケージを (ドキュメントに記載されているように) 「cheapo 空間データベース」として使用することをお勧めします。考え方は主に次のようになります。
#!/usr/bin/env python
from itertools import product
from random import uniform, sample, seed
from rtree import index
from shapely.geometry import Point, Polygon, box, shape
from shapely.affinity import translate
seed(666)
#generate random polygons, in your case, the polygons are stored
#in geo_data['features']
P = Polygon([(0, 0), (0.5, 0), (0.5, 0.5), (0, 0.5), (0, 0)])
polygons = []
for dx, dy in product(range(0, 100), range(0, 100)):
polygons.append(translate(P, dx, dy))
#construct the spatial index and insert bounding boxes of all polygons
idx = index.Index()
for pid, P in enumerate(polygons):
idx.insert(pid, P.bounds)
delta = 0.5
for i in range(0, 1000):
#generate random points
x, y = uniform(0, 10), uniform(0, 10)
pnt = Point(x, y)
#create a region around the point of interest
bounds = (x-delta, y-delta, x+delta, y+delta)
#also possible, but much slower
#bounds = pnt.buffer(delta).bounds
#the index tells us which polygons are worth checking, i.e.,
#the bounding box of which intersects with the region constructed in previous step
for candidate in idx.intersection(bounds):
P = polygons[candidate]
#test only these candidates
if P.contains(pnt):
print(pnt, P)
この記事はインターネットから収集されたものであり、転載の際にはソースを示してください。
侵害の場合は、連絡してください[email protected]
コメントを追加