Finden von Nachbarn erster Ordnung mit Shapefile-Polygonen

Question

Ich suche eine effiziente Möglichkeit, die Nachbarn 1. Ordnung eines gegebenen Polygons zu finden. Meine Daten sind im Format shapefile.

Meine erste Idee war, die x- und y-Koordinaten der Polygone zu berechnen, um die Zentroide des Nachbarn zu finden.

import pysal 
from pysal.common import * 
import pysal.weights 
import numpy as np 
from scipy import sparse,float32 
import scipy.spatial 
import os, gc, operator 


def get_points_array_from_shapefile(inFile): 
    """ 
    Gets a data array of x and y coordinates from a given shape file 

    Parameters 
    ---------- 
    shapefile: string name of a shape file including suffix 

    Returns 
    ------- 
    points: array (n,2) a data array of x and y coordinates 

    Notes 
    ----- 
    If the given shape file includes polygons, 
    this function returns x and y coordinates of the polygons' centroids 

    Examples 
    -------- 
    Point shapefile 
    >>> from pysal.weights.util import get_points_array_from_shapefile 
    >>> xy = get_points_array_from_shapefile('../examples/juvenile.shp') 
    >>> xy[:3] 
    array([[ 94., 93.], 
      [ 80., 95.], 
      [ 79., 90.]]) 

    Polygon shapefile 
    >>> xy = get_points_array_from_shapefile('../examples/columbus.shp') 
    >>> xy[:3] 
    array([[ 8.82721847, 14.36907602], 
      [ 8.33265837, 14.03162401], 
      [ 9.01226541, 13.81971908]]) 

    (source: https://code.google.com/p/pysal/source/browse/trunk/pysal/weights/util.py?r=1013) 

    """ 
    f = pysal.open(inFile) 
    shapes = f.read() 
    if f.type.__name__ == 'Polygon': 
     data = np.array([shape.centroid for shape in shapes]) 
    elif f.type.__name__ == 'Point': 
     data = np.array([shape for shape in shapes]) 
    f.close() 
    return data 


inFile = "../examples/myshapefile.shp" 
my_centr = get_points_array_from_shapefile(inFile) 

Dieser Ansatz könnte für ein reguläres Raster gelten, aber in meinem Fall muss ich eine "allgemeinere" Lösung finden. Die Abbildung zeigt das Problem. Betrachte das gelbe Polygon als Schiedsrichter. Die Polygone des Nachbarn sind die grauen Polygone. Bei Verwendung des Schwerpunkts-Nachbarn-Ansatzes wird das klare blaue Polygon als Nachbar betrachtet, es hat jedoch keine gemeinsame Seite mit dem gelben Polygon.

Eine neuere Lösung aus Efficiently finding the 1st order neighbors of 200k polygons modifiziert, um die folgenden sein können:

from collections import defaultdict 
inFile = 'C:\\MultiShapefile.shp' 

shp = osgeo.ogr.Open(inFile) 
layer = shp.GetLayer() 
BlockGroupVertexDictionary = dict() 
for index in xrange(layer.GetFeatureCount()): 
    feature = layer.GetFeature(index) 
    FID = str(feature.GetFID()) 
    geometry = feature.GetGeometryRef() 
    pts = geometry.GetGeometryRef(0) 
    # delete last points because is the first (see shapefile polygon topology) 
    for p in xrange(pts.GetPointCount()-1): 
     PointText = str(pts.GetX(p))+str(pts.GetY(p)) 
     # If coordinate is already in dictionary, append this BG's ID 
     if PointText in BlockGroupVertexDictionary: 
      BlockGroupVertexDictionary[PointText].append(FID) 
     # If coordinate is not already in dictionary, create new list with this BG's ID 
     else: 
      BlockGroupVertexDictionary[PointText] = [FID] 

Mit dieser Lösung haben ein Wörterbuch mit dem Scheitelpunkt als die Tasten-Koordinaten und eine Liste der Blockgruppe IDs, die einen Scheitelpunkt an, dass haben Koordinaten als der Wert.

>>> BlockGroupVertexDictionary 
{'558324.3057036361423.57178': ['18'], 
'558327.4401686361422.40755': ['18', '19'], 
'558347.5890836361887.12271': ['1'], 
'558362.8645026361662.38757': ['17', '18'], 
'558378.7836876361760.98381': ['14', '17'], 
'558389.9225016361829.97259': ['14'], 
'558390.1235856361830.41498': ['1', '14'], 
'558390.1870856361652.96599': ['17', '18', '19'], 
'558391.32786361398.67786': ['19', '20'], 
'558400.5058556361853.25597': ['1'], 
'558417.6037156361748.57558': ['14', '15', '17', '19'], 
'558425.0594576362017.45522': ['1', '3'], 
'558438.2518686361813.61726': ['14', '15'], 
'558453.8892486362065.9571': ['3', '5'], 
'558453.9626046361375.4135': ['20', '21'], 
'558464.7845966361733.49493': ['15', '16'], 
'558474.6171066362100.82867': ['4', '5'], 
'558476.3606496361467.63697': ['21'], 
'558476.3607186361467.63708': ['26'], 
'558483.1668826361727.61931': ['19', '20'], 
'558485.4911846361797.12981': ['15', '16'], 
'558520.6376956361649.94611': ['25', '26'], 
'558525.9186066361981.57914': ['1', '3'], 
'558527.5061096362189.80664': ['4'], 
'558529.0036896361347.5411': ['21'], 
'558529.0037236361347.54108': ['26'], 
'558529.8873646362083.17935': ['4', '5'], 
'558533.062376362006.9792': ['1', '3'], 
'558535.4436256361710.90985': ['9', '16', '20'], 
'558535.4437266361710.90991': ['25'], 
'558548.7071816361705.956': ['9', '10'], 
'558550.2603156361432.56769': ['26'], 
'558550.2603226361432.56763': ['21'], 
'558559.5872216361771.26884': ['9', '16'], 
'558560.3288756362178.39003': ['4', '5'], 
'558568.7811926361768.05997': ['1', '9', '10'], 
'558572.749956362041.11051': ['3', '5'], 
'558573.5437016362012.53546': ['1', '3'], 
'558575.3048386362048.77518': ['2', '3'], 
'558576.189546362172.87328': ['5'], 
'558577.1149386361695.34587': ['7', '10'], 
'558579.0999636362020.47297': ['1', '3'], 
'558581.6312396362025.36096': ['0', '1'], 
'558586.7728956362035.28967': ['0', '3'], 
'558589.8015336362043.7987': ['2', '3'], 
'558601.3250076361686.30355': ['7'], 
'558601.3250736361686.30353': ['25'], 
'558613.7793476362164.19871': ['2', '5'], 
'558616.4062876361634.7097': ['7'], 
'558616.4063116361634.70972': ['25'], 
'558618.129066361634.29952': ['7', '11', '22'], 
'558618.1290896361634.2995': ['25'], 
'558626.9644156361875.47515': ['10', '11'], 
'558631.2229836362160.17325': ['2'], 
'558632.0261236361600.77448': ['25', '26'], 
'558639.495586361898.60961': ['11', '13'], 
'558650.4935686361918.91358': ['12', '13'], 
'558659.2473416361624.50945': ['8', '11', '22', '24'], 
'558664.5218136361857.94836': ['7', '10'], 
'558666.4126376361622.80343': ['8', '24'], 
'558675.1439056361912.52276': ['12', '13'], 
'558686.3385396361985.08892': ['0', '1'], 
.................. 
................. 
'558739.4377836361931.57279': ['11', '13'], 
'558746.8758486361973.84475': ['11', '13'], 
'558751.3440576361902.20399': ['6', '11'], 
'558768.8067026361258.4715': ['26'], 
'558779.9170276361961.16408': ['6', '11'], 
'558785.7399596361571.47416': ['22', '24'], 
'558791.5596546361882.09619': ['8', '11'], 
'558800.2351726361877.75843': ['6', '8'], 
'558802.7700816361332.39227': ['26'], 
'558802.770176361332.39218': ['22'], 
'558804.7899976361336.78827': ['22'], 
'558812.9707376361565.14513': ['23', '24'], 
'558833.2667696361940.68932': ['6', '24'], 
'558921.2068976361539.98868': ['22', '23'], 
'558978.3570116361885.00604': ['23', '24'], 
'559022.80716361982.3729': ['23'], 
'559096.8905816361239.42141': ['22'], 
'559130.7573166361935.80614': ['23'], 
'559160.3907086361434.15513': ['22']} 

Answer 1

bin ich mit den spezifischen Datenformate verwendet werden, nicht vertraut aber egal, denken würde, die folgende Idee arbeiten.

In Python können Sie Sätze aus Tupeln von Zahlen machen, d. H. (x,y) und (x1,y1,x2,y2), so dass es möglich sein sollte, einen Satz zu machen, der alle Punkte oder Kanten in einem gegebenen Polygon repräsentiert. Danach können Sie sehr schnelle Kreuzungsoperationen verwenden, um alle Nachbarn erster Ordnung zu finden.

Sie könnten in der Lage sein, den Prozess zu beschleunigen, eine Art trivial Ablehnung Test mit der Weiterverarbeitung von Polygonen zu vermeiden, die möglicherweise nicht ein Nachbar sein könnten - vielleicht Ihren Polygone Centroide Idee verwenden.

Macht diese Erklärung Sinn?

Answer 2

Nur wenn dies noch eine offene Frage für das OP ist oder jemand anderes stolpert hier.

import pysal as ps 
w = ps.queen_from_shapefile('shapefile.shp') 

http://pysal.readthedocs.io/en/latest/users/tutorials/weights.html#pysal-spatial-weight-types