寻找一种快速的方法来使用 Shapely 找到一个点所属的多边形 [英] Looking for a fast way to find the polygon a point belongs to using Shapely

问题描述

我有一组 ~36,000 个多边形，它们代表了该国的一个分区(~县).我的python脚本接收了很多点:pointId、经度、纬度.

I have a set of ~36,000 polygons which represent a partition (~counties) of the country. My python script receives a lot of points: pointId, longitude, latitude.

对于每个点，我想发回 pointId、polygonId.对于每个点，循环到所有多边形并使用 myPoint.within(myPolygon) 效率很低.

For each point, I want to send back pointId, polygonId. For each point, looping into all the polygons and use myPoint.within(myPolygon) is quite inefficient.

我认为 shapely 库提供了一种更好的方法来准备多边形，以便为点查找多边形成为树路径(国家、地区、子地区，...)

I suppose the shapely library offers a better way to prepare the polygon so that finding the polygon for a point becomes a tree path (country, region, sub region, ...)

这是我目前的代码:

import sys
import os
import json
import time
import string
import uuid

py_id = str(uuid.uuid4())

sys.stderr.write(py_id + '
')
sys.stderr.write('point_in_polygon.py V131130a.
')
sys.stderr.flush()

from shapely.geometry import Point
from shapely.geometry import Polygon
import sys
import json
import string
import uuid
import time

jsonpath='.cantons.json'
jsonfile = json.loads(open(jsonpath).read())

def find(id, obj):
    results = []

    def _find_values(id, obj):
        try:
            for key, value in obj.iteritems():
                if key == id:
                    results.append(value)
                elif not isinstance(value, basestring):
                    _find_values(id, value)
        except AttributeError:
            pass

        try:
            for item in obj:
                if not isinstance(item, basestring):
                    _find_values(id, item)
        except TypeError:
            pass

    if not isinstance(obj, basestring):
        _find_values(id, obj)
    return results


#1-Load polygons from json  
r=find('rings',jsonfile)
len_r = len(r)

#2-Load attributes from json
a=find('attributes',jsonfile)

def insidePolygon(point,json):
        x=0
    while x < len_r :
            y=0
            while y < len(r[x]) :
        p=Polygon(r[x][y])
        if(point.within(p)):
            return a[y]['OBJECTID'],a[y]['NAME_5']
        y=y+1
        x=x+1
    return None,None


while True:
    line = sys.stdin.readline()
    if not line:
        break
    try:
        args, tobedropped = string.split(line, "
", 2)

        #input: vehicleId, longitude, latitude
        vehicleId, longitude, latitude = string.split(args, "	")

        point = Point(float(longitude), float(latitude))
        cantonId,cantonName = insidePolygon(point,r)

        #output: vehicleId, cantonId, cantonName
        # vehicleId will be 0 if not found
        # vehicleId will be < 0 in case of an exception
        if (cantonId == None):
            print "	".join(["0", "", ""])
        else:
            print "	".join([str(vehicleId), str(cantonId), str(cantonName)])

    except ValueError:
        print "	".join(["-1", "", ""])
        sys.stderr.write(py_id + '
')
        sys.stderr.write('ValueError in Python script
')
        sys.stderr.write(line)
        sys.stderr.flush()
    except:
        sys.stderr.write(py_id + '
')
        sys.stderr.write('Exception in Python script
')
        sys.stderr.write(str(sys.exc_info()[0]) + '
')
        sys.stderr.flush()
        print "	".join(["-2", "", ""])

推荐答案

使用 Rtree (examples) 为 R 树索引 以:(1)索引 36k 多边形的边界(在读取 jsonfile 后执行此操作)，然后(2)非常快速地找到相交的边界框每个多边形到您的兴趣点.然后，(3)对于来自 Rtree 的相交边界框，使用 shapely 来使用，例如point.within(p) 进行实际的多边形点分析.使用这种技术，您应该会看到性能的巨大飞跃.

Use Rtree (examples) as R-tree index to: (1) index the bounds of the 36k polygons (do this just after reading jsonfile), then (2) very quickly find the intersecting bounding boxes of each polygon to your point of interest. Then, (3) for the intersecting bounding boxes from Rtree, use shapely to use, e.g. point.within(p) to do the actual point-in-polygon analysis. You should see a massive leap of performance with this technique.

这篇关于寻找一种快速的方法来使用 Shapely 找到一个点所属的多边形的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！