Я попробовал несколько примеров кода с использованием библиотек, таких как shapefile, fiona и ogr, чтобы попытаться проверить, попадает ли точка (x, y) в границы мультиполигона, созданного с помощью ArcMap (и, следовательно, в формате shapefile). Однако ни один из примеров не работает хорошо с мультиполигонами, хотя они хорошо работают с обычными, однополигональными шейп-файлами. Ниже приведены некоторые фрагменты, которые я попробовал:
# 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'
Первый пример отлично работает с одним полигоном за раз, но когда я ввожу точку в пределах одной из фигур в моем многогранном шейп-файле, он возвращает «out», даже если он попадает в одну из множества частей.
Во втором примере я получаю сообщение об ошибке «объект типа« Геометрия »не имеет len ()», который, как я полагаю, связан с тем, что поле геометрии не может быть прочитано как обычный индексированный список / массив.
Я также попытался заменить фактическую точку в коде многоугольника, как предложено здесь, чтобы убедиться, что это не та часть кода, которая не работает. И хотя примеры этой ссылки хорошо работают с простыми шейп-файлами полигонов, я не могу заставить мой сложный мультиполигон правильно тестировать.
Поэтому я не могу придумать какой-либо другой способ проверить, попадает ли точка в шейп-файл многоугольника через python ... Может быть, есть другие библиотеки, которые мне не хватает?
источник
polygon = Polygon(geometry)
какой-то цикл try, где он переключается,polygon = MultiPolygon(geometry)
если возникает эта ошибка.Ответы:
Шейп-файлы не имеют типа MultiPolygon (type = Polygon), но они все равно поддерживают их (все кольца хранятся в одном элементе = список полигонов, посмотрите на Преобразование огромного мультиполигона в полигоны )
Проблема
Если я открою шейп-файл MultiPolygon, геометрия будет «Полигон»
Решение 1 с Фионой
Fiona интерпретирует эту функцию как MultiPolygon, и вы можете применить решение, представленное в разделе « Более эффективное пространственное соединение» в Python, без QGIS, ArcGIS, PostGIS и т. Д. (1)
Решение 2 с pyshp (shapefile) и протоколом geo_interface (как GeoJSON)
Это дополнение к ответу xulnik.
Решение 3 с ogr и протоколом geo_interface ( приложения Python Geo_interface )
Решение 4 с GeoPandas, как в более эффективном пространственном соединении в Python без QGIS, ArcGIS, PostGIS и т. Д. (2)
Точки 1,3,5,6 попадают в границы мультиполигона
источник
multi = shape(polys.shapeRecords()[0].shape.__geo_interface__)
Решение 2? Я не могу получить вызов метода shape () изshapefile.py
. Я даже пыталсяshapefile.Shape()
; есть класс для него, но он не работает.within()
метод?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'
Проблема в вашем первом примере заключается в следующем цикле:
Добавляет только последние характерные точки. Я опробовал свой подход с этим шейп-файлом:
Я изменил ваш код:
Вышеприведенный код был запущен на консоли Python QGIS, и в результате было получено:
Он работает отлично, и теперь вы можете проверить, попадает ли точка (x, y) в границы каждого объекта.
источник
Если вы пытаетесь проверить широту, долготу точки внутри многоугольника, убедитесь, что у вас есть точечный объект, созданный с помощью следующего:
Точка принимает долготу, затем широту в аргументе. Сначала не широта. Вы можете вызвать
polygon_object.within
функцию, чтобы проверить, находится ли точка внутри фигуры.источник