Поиск соседей 1-го порядка с использованием полигонов шейп-файла

Я ищу эффективный способ найти соседей 1-го порядка данного многоугольника. Мои данные находятся в формате shapefile.

Моей первой идеей было вычислить координаты x и y центроидов полигонов, чтобы найти центроиды соседей.

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)

Этот подход может быть применим для обычной сетки, но в моем случае мне нужно найти «более общее» решение. На рисунке показана проблема. Учтите, что у желтого многоугольника есть судья. Соседние полигоны — серые полигоны. При использовании подхода центроиды-соседи чистый синий многоугольник считается соседним, но у него нет общей стороны с желтым многоугольником.

Недавнее решение, измененное из эффективного поиска соседей 1-го порядка для 200 000 полигонов могут быть следующими:

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]

С этим решением у меня есть словарь с координатами вершин в качестве ключей и список идентификаторов группы блоков, у которых есть вершина с этой координатой в качестве значения.

>>> 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']}

введите здесь описание изображения


person Gianni Spear    schedule 21.11.2012    source источник
comment
Каково точное определение соседа 1-го порядка на нерегулярной сетке? Означает ли это просто наличие общего преимущества?   -  person martineau    schedule 21.11.2012
comment
@martineau, соседи 1-го порядка - это все многоугольники с общей границей (= вершина в случае шейп-файла) с многоугольником-i   -  person Gianni Spear    schedule 21.11.2012
comment
@мартино. Я нахожу эту ссылку в гугле. gis.stackexchange.com/questions/17457/ выглядит хорошей стратегией, но я хочу работать вне модуля Arcmap   -  person Gianni Spear    schedule 21.11.2012
comment
Посмотрите на стройные. toblerity.github.com/shapely/manual.html   -  person Joe Kington    schedule 21.11.2012
comment
@ Джанни, не уверен, что ты когда-нибудь получил ответ, но если ты хочешь сделать это в R, есть довольно простая функция, называемая poly2nb, как часть пакета spdep: inside-r.org/packages/cran/spdep/docs/poly2nb. Если бы вы получили ответ на Python, мне было бы очень интересно его узнать.   -  person maxliving    schedule 08.04.2014
comment
Если вы хотите загрузить свои пространственные веса в библиотеку pysal, существует метод pysal.shimbel, который для заданных весов возвращает отношения первого, второго, третьего порядка (и т. д.) для каждого наблюдения.   -  person Sam Comber    schedule 01.05.2018


Ответы (2)


На всякий случай это все еще открытый вопрос для ОП или кто-то еще спотыкается здесь.

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

person Jzl5325    schedule 16.01.2015

Я не знаком с конкретными используемыми форматами данных, но, тем не менее, думаю, что следующая идея сработает.

В Python вы можете создавать наборы из кортежей чисел, то есть (x,y) и (x1,y1,x2,y2), поэтому должна быть возможность создать набор, представляющий все точки или ребра в заданном многоугольнике. После этого вы сможете использовать очень быстрые операции пересечения множества, чтобы найти всех соседей 1-го порядка.

Возможно, вы сможете ускорить процесс, используя какой-нибудь тривиальный тест отбраковки, чтобы избежать дальнейшей обработки полигонов, которые не могут быть соседними — возможно, используя идею центроидов ваших полигонов.

Имеет ли это объяснение смысл?

person martineau    schedule 21.11.2012
comment
так себе, извините. Я пытаюсь написать новый подход, чтобы хранить в словаре все многоугольники с одной вершиной за многоугольником и добавлять идентификатор каждого многоугольника. - person Gianni Spear; 21.11.2012
comment
@Gianni: Неясно, какое отношение это имеет к поиску соседей 1-го порядка данного многоугольника - вам придется сравнивать вершины или векторы ребер на каком-то уровне, чтобы сделать это. - person martineau; 21.11.2012
comment
взгляните на мое новое решение. Из словаря мне нужно понять, например, многоугольник «18» многоугольники с общими вершинами - person Gianni Spear; 21.11.2012
comment
Множество чем-то похоже на словарь (на самом деле когда-то они были реализованы с их использованием), но они были оптимизированы для операций над множествами, таких как объединение, пересечение, разность и т. д. Несмотря на это, я предлагаю одну вещь: это было бы более эффективно с обоими набор или словарь для использования ключей, подобных этому (pts.GetX(p), pts.GetY(p)), который является кортежем вместо str(pts.GetX(p))+str(pts.GetY(p)), используемого в только что добавленном коде решения. - person martineau; 21.11.2012