Próbowałem kilku przykładów kodu wykorzystującego biblioteki takie jak shapefile, fiona i ogr, aby spróbować sprawdzić, czy punkt (x, y) mieści się w granicach wieloboku utworzonego za pomocą ArcMap (a zatem w formacie shapefile). Jednak żaden z przykładów nie działa dobrze z multipoligonami, chociaż dobrze sobie radzą z normalnymi, pojedynczymi wielokątami kształtowymi. Niektóre fragmenty, które wypróbowałem, są poniżej:
# First example using shapefile and shapely:
from shapely.geometry import Polygon, Point, MultiPolygon
import shapefile
polygon = shapefile.Reader('shapefile.shp')
polygon = polygon.shapes()
shpfilePoints = []
for shape in polygon:
shpfilePoints = shape.points
polygon = shpfilePoints
poly = Polygon(poly)
point = Point(x, y)
# point in polygon test
if polygon.contains(point):
print 'inside'
else:
print 'OUT'
# Second example using ogr and shapely:
from shapely.geometry import Polygon, Point, MultiPolygon
from osgeo import ogr, gdal
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open("shapefile.shp", 0)
layer = dataset.GetLayer()
for index in xrange(layer.GetFeatureCount()):
feature = layer.GetFeature(index)
geometry = feature.GetGeometryRef()
polygon = Polygon(geometry)
print 'polygon points =', polygon # this prints 'multipoint' + all the points fine
point = Point(x, y)
# point in polygon test
if polygon.contains(point):
print 'inside'
else:
print 'OUT'
Pierwszy przykład działa dobrze z pojedynczym wielokątem na raz, ale kiedy wprowadzam punkt w jednym z kształtów w moim wielobokowym pliku kształtu, zwraca on „out”, nawet jeśli mieści się w jednej z wielu części.
W drugim przykładzie pojawia się błąd „obiekt typu„ Geometria ”nie ma len ()”, co, jak zakładam, jest spowodowane tym, że pola geometrii nie można odczytać jako normalnej, indeksowanej listy / tablicy.
Dodatkowo próbowałem zastąpić rzeczywisty punkt w kodzie wielokąta, jak sugerowano tutaj, aby upewnić się, że nie jest to ta część kodu, która nie działa. I chociaż przykłady tego linku działają dobrze z prostymi kształtami wielokątów, nie mogę poprawnie przetestować mojego złożonego multipoligonu.
Nie mogę więc wymyślić żadnego innego sposobu na sprawdzenie, czy punkt mieści się w pliku kształtu wielokąta za pośrednictwem pytona ... Może istnieją inne biblioteki, których mi brakuje?
źródło
polygon = Polygon(geometry)
jakąś pętlą try, w której przełącza się,polygon = MultiPolygon(geometry)
jeśli wystąpi ten błąd.Odpowiedzi:
Pliki kształtów nie mają typu MultiPolygon (type = Polygon), ale i tak je obsługują (wszystkie pierścienie są przechowywane w jednej funkcji = lista wielokątów, spójrz na Konwertowanie ogromnego wieloboku na wieloboki )
Problem
Jeśli otworzę plik kształtu MultiPolygon, geometria to „Wielokąt”
Rozwiązanie 1 z Fioną
Fiona interpretuje tę funkcję jako MultiPolygon i możesz zastosować rozwiązanie przedstawione w Bardziej wydajnym złączeniu przestrzennym w Pythonie bez QGIS, ArcGIS, PostGIS itp. (1)
Rozwiązanie 2 z pyshp (shapefile) i protokołem geo_interface (podobnym do GeoJSON)
To jest dodatek do odpowiedzi xulnika.
Rozwiązanie 3 z ogrodem i protokołem geo_interface ( aplikacje Python Geo_interface )
Rozwiązanie 4 z GeoPandas jak w bardziej wydajnym sprzężeniu przestrzennym w Pythonie bez QGIS, ArcGIS, PostGIS itp. (2)
Punkty 1,3,5,6 mieszczą się w granicach MultiPolygon
źródło
multi = shape(polys.shapeRecords()[0].shape.__geo_interface__)
w Solution 2? Nie mogę uzyskać wywołania metody shape () zshapefile.py
. Próbowałem nawetshapefile.Shape()
; jest na to klasa, ale to nie działa.within()
metodę?from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon
)File "C:\WinPython\python-3.6.5.amd64\lib\site-packages\geopandas\tools\sjoin.py", line 43, in sjoin if left_df.crs != right_df.crs:
AttributeError: 'MultiPolygon' object has no attribute 'crs'
Problem w twoim pierwszym przykładzie dotyczy tej pętli:
Dołącza tylko ostatnie punkty funkcji. Wypróbowałem swoje podejście z tym plikiem kształtu:
Zmodyfikowałem twój kod do:
Powyższy kod został uruchomiony w konsoli Pythona w QGIS, a wynik był następujący:
Działa idealnie i teraz możesz sprawdzić, czy punkt (x, y) mieści się w granicach każdej cechy.
źródło
Jeśli próbujesz sprawdzić szerokość i długość geograficzną w wielokącie, upewnij się, że masz obiekt punktowy utworzony przez:
Punkt ma długość geograficzną, a następnie szerokość geograficzną w argumencie. Najpierw szerokość geograficzna. Możesz wywołać
polygon_object.within
funkcję, aby sprawdzić, czy punkt znajduje się w kształcie.źródło