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

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?

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”

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

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:

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:

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

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.