Znajdowanie sąsiadów pierwszego rzędu za pomocą kształtów wielokątów

Question

Szukam skutecznego sposobu na znalezienie sąsiadów pierwszego rzędu danego wielokąta. Moje dane są w formacie shapefile.

Mój pierwszy pomysł polegał na obliczeniu współrzędnych xiy współrzędnych wielokątów w celu znalezienia centroidów sąsiada.

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) 

Takie podejście może być ważne dla regularnej siatki, ale w moim przypadku muszę znaleźć "bardziej ogólne" rozwiązanie. Rysunek pokazuje problem. Zastanów się, czy żółty wielokąt ma sędziego. Poligony sąsiada są szarymi wielokątami. Używając centroidów-sąsiadów podejście, przezroczysty niebieski wielokąt jest uważany za sąsiada, ale nie ma strony wspólnej z żółtym wielokątem.

Ostatnie rozwiązanie modyfikowany Efficiently finding the 1st order neighbors of 200k polygons mogą być następujące:

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] 

Dzięki takiemu rozwiązaniu, mam słownik z wierzchołkiem jako zakres kluczy i listę identyfikatorów grupy bloków, które mają wierzchołek na to współrzędna jako wartość.

>>> 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

nie jestem zaznajomiony z używane specyficzne formaty danych, ale niezależnie od tego, myślę następujące pomysł będzie działać.

W języku Python można tworzyć zestawy z krotek liczb, tj. (x,y) i (x1,y1,x2,y2), dlatego powinno być możliwe utworzenie zestawu reprezentującego wszystkie punkty lub krawędzie w danym wielokącie. Następnie będziesz mógł użyć bardzo szybkiego zestawu operacji przecięcia, aby znaleźć wszystkich sąsiadów pierwszego rzędu.

Być może będziesz w stanie przyspieszyć proces, stosując jakiś trywialny test odrzucania, aby uniknąć dalszego przetwarzania wielokątów, które nie mogłyby być sąsiadem - być może używając idei centroidów wielokątów.

Czy to wyjaśnienie ma sens?

Answer 2

To tylko otwarte pytanie do OP, albo ktoś tu się potknie.

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