Sprawdź, czy punkt mieści się w wieloboku z Pythonem

13

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?

spartmar
źródło
twój drugi przykład wygląda na przymuszanie wieloboku do wielokąta? Może to być tylko sprawdzanie punktu względem pierwszej części wieloboku. Spróbuj przenieść punkt do różnych części i sprawdź, czy sprawdzenie kiedykolwiek się powiedzie.
obrl_soil
@obrl_soil Dziękujemy za sugestię. Jednak drugi przykład nigdy nie działa z powodu komunikatu o błędzie, który opisałem powyżej (obiekt typu „Geometria” nie ma len ()) „bez względu na to, czy próbuję MultiPolygon (geometria) czy po prostu Polygon (geometria). Próbowałem wielu punktów w pierwszy przykład i tylko ci z głównego dzieła w zakresie wielokątów. Mam nadzieję, że to wyjaśnienie pomoże
spartmar
Tak, myślę, że musisz zastąpić polygon = Polygon(geometry)jakąś pętlą try, w której przełącza się, polygon = MultiPolygon(geometry)jeśli wystąpi ten błąd.
obrl_soil
Problem w twoim pierwszym przykładzie dotyczy pierwszej pętli.
xunilk

Odpowiedzi:

24

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

wprowadź opis zdjęcia tutaj

Jeśli otworzę plik kształtu MultiPolygon, geometria to „Wielokąt”

multipolys = fiona.open("multipol.shp")
multipolys.schema
{'geometry': 'Polygon', 'properties': OrderedDict([(u'id', 'int:10')])}
len(multipolys)
1

Rozwiązanie 1 z Fioną

import fiona
from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon
multipol = fiona.open("multipol.shp")
multi= multipol.next() # only one feature in the shapefile
print multi
{'geometry': {'type': 'MultiPolygon', 'coordinates': [[[(-0.5275288092189501, 0.5569782330345711), (-0.11779769526248396, 0.29065300896286816), (-0.25608194622279135, 0.01920614596670933), (-0.709346991037132, -0.08834827144686286), (-0.8629961587708066, 0.18309859154929575), (-0.734955185659411, 0.39820742637644047), (-0.5275288092189501, 0.5569782330345711)]], [[(0.19974391805377723, 0.060179257362355965), (0.5480153649167734, 0.1293213828425096), (0.729833546734955, 0.03969270166453265), (0.8143405889884763, -0.13956466069142115), (0.701664532650448, -0.38540332906530095), (0.4763124199743918, -0.5006402048655569), (0.26888604353393086, -0.4238156209987196), (0.18950064020486557, -0.2291933418693981), (0.19974391805377723, 0.060179257362355965)]], [[(-0.3764404609475033, -0.295774647887324), (-0.11523687580025621, -0.3597951344430217), (-0.033290653008962945, -0.5800256081946222), (-0.11523687580025621, -0.7413572343149808), (-0.3072983354673495, -0.8591549295774648), (-0.58898847631242, -0.6927016645326505), (-0.6555697823303457, -0.4750320102432779), (-0.3764404609475033, -0.295774647887324)]]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', 1)])}

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)

points= ([pt for pt  in fiona.open("points.shp")])
for i, pt in enumerate(points):
    point = shape(pt['geometry'])
    if point.within(shape(multi['geometry'])):
         print i, shape(points[i]['geometry'])
1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.4993597951344431 -0.06017925736235585)
5 POINT (-0.3764404609475033 -0.4750320102432779)
6 POINT (-0.3098591549295775 -0.6312419974391805)

Rozwiązanie 2 z pyshp (shapefile) i protokołem geo_interface (podobnym do GeoJSON)

To jest dodatek do odpowiedzi xulnika.

import shapefile
pts = shapefile.Reader("points.shp")
polys = shapefile.Reader("multipol.shp")
points = [pt.shape.__geo_interface__ for pt in pts.shapeRecords()]
multi = shape(polys.shapeRecords()[0].shape.__geo_interface__) # 1 polygon
print multi
MULTIPOLYGON (((-0.5275288092189501 0.5569782330345711, -0.117797695262484 0.2906530089628682, -0.2560819462227913 0.01920614596670933, -0.7093469910371319 -0.08834827144686286, -0.8629961587708066 0.1830985915492958, -0.734955185659411 0.3982074263764405, -0.5275288092189501 0.5569782330345711)), ((0.1997439180537772 0.06017925736235596, 0.5480153649167734 0.1293213828425096, 0.729833546734955 0.03969270166453265, 0.8143405889884763 -0.1395646606914211, 0.701664532650448 -0.3854033290653009, 0.4763124199743918 -0.5006402048655569, 0.2688860435339309 -0.4238156209987196, 0.1895006402048656 -0.2291933418693981, 0.1997439180537772 0.06017925736235596)), ((-0.3764404609475033 -0.295774647887324, -0.1152368758002562 -0.3597951344430217, -0.03329065300896294 -0.5800256081946222, -0.1152368758002562 -0.7413572343149808, -0.3072983354673495 -0.8591549295774648, -0.58898847631242 -0.6927016645326505, -0.6555697823303457 -0.4750320102432779, -0.3764404609475033 -0.295774647887324)))
for i, pt in enumerate(points):
    point = shape(pt)
    if point.within(multi): 
        print i, shape(points[i])
1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.4993597951344431 -0.06017925736235585)
5 POINT (-0.3764404609475033 -0.4750320102432779)
6 POINT (-0.3098591549295775 -0.6312419974391805)

Rozwiązanie 3 z ogrodem i protokołem geo_interface ( aplikacje Python Geo_interface )

from osgeo import ogr
import json
def records(file):  
    # generator 
    reader = ogr.Open(file)
    layer = reader.GetLayer(0)
    for i in range(layer.GetFeatureCount()):
        feature = layer.GetFeature(i)
        yield json.loads(feature.ExportToJson())

points  = [pt for pt in records("point_multi_contains.shp")]
multipol = records("multipol.shp")
multi = multipol.next() # 1 feature
for i, pt in enumerate(points):
     point = shape(pt['geometry'])
     if point.within(shape(multi['geometry'])):
          print i, shape(points[i]['geometry'])

1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.499359795134443 -0.060179257362356)
5 POINT (-0.376440460947503 -0.475032010243278)
6 POINT (-0.309859154929577 -0.631241997439181)

Rozwiązanie 4 z GeoPandas jak w bardziej wydajnym sprzężeniu przestrzennym w Pythonie bez QGIS, ArcGIS, PostGIS itp. (2)

import geopandas
point = geopandas.GeoDataFrame.from_file('points.shp') 
poly  = geopandas.GeoDataFrame.from_file('multipol.shp')
from geopandas.tools import sjoin
pointInPolys = sjoin(point, poly, how='left')
grouped = pointInPolys.groupby('index_right')
list(grouped)
[(0.0,      geometry                               id_left  index_right id_right  

1  POINT (-0.58898847631242 0.17797695262484)       None      0.0        1.0 
3  POINT (0.4993597951344431 -0.06017925736235585)  None      0.0        1.0
5  POINT (-0.3764404609475033 -0.4750320102432779)  None      0.0        1.0 
6  POINT (-0.3098591549295775 -0.6312419974391805)  None      0.0        1.0 ]
print grouped.groups
{0.0: [1, 3, 5, 6]} 

Punkty 1,3,5,6 mieszczą się w granicach MultiPolygon

gen
źródło
Nieco stary wątek tutaj, ale jak wywołać multi = shape(polys.shapeRecords()[0].shape.__geo_interface__)w Solution 2? Nie mogę uzyskać wywołania metody shape () z shapefile.py. Próbowałem nawet shapefile.Shape(); jest na to klasa, ale to nie działa.
pstatix
Ponadto, skąd bierzesz tę within()metodę?
pstatix
1
z Shapely ( from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon)
gen
Ten błąd 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'
Aaron Bramson,
6

Problem w twoim pierwszym przykładzie dotyczy tej pętli:

...
shpfilePoints = []
for shape in polygon:
    shpfilePoints = shape.points
...

Dołącza tylko ostatnie punkty funkcji. Wypróbowałem swoje podejście z tym plikiem kształtu:

wprowadź opis zdjęcia tutaj

Zmodyfikowałem twój kod do:

from shapely.geometry import Polygon, Point, MultiPolygon
import shapefile 

path = '/home/zeito/pyqgis_data/polygon8.shp'

polygon = shapefile.Reader(path) 

polygon = polygon.shapes() 

shpfilePoints = [ shape.points for shape in polygon ]

print shpfilePoints

polygons = shpfilePoints

for polygon in polygons:
    poly = Polygon(polygon)
    print poly

Powyższy kod został uruchomiony w konsoli Pythona w QGIS, a wynik był następujący:

wprowadź opis zdjęcia tutaj

Działa idealnie i teraz możesz sprawdzić, czy punkt (x, y) mieści się w granicach każdej cechy.

Xunilk
źródło
0

Jeśli próbujesz sprawdzić szerokość i długość geograficzną w wielokącie, upewnij się, że masz obiekt punktowy utworzony przez:

from shapely.geometry.point import Point
Point(LONGITUDE, LATITUDE)
..
poly.within(point) # Returns true if the point within the 

Punkt ma długość geograficzną, a następnie szerokość geograficzną w argumencie. Najpierw szerokość geograficzna. Możesz wywołać polygon_object.withinfunkcję, aby sprawdzić, czy punkt znajduje się w kształcie.

Ahmed
źródło