Как я могу проверить, находится ли геоточка в области данного шейп-файла?
Мне удалось загрузить шейп-файл в python, но я не могу продвинуться дальше.
Как я могу проверить, находится ли геоточка в области данного шейп-файла?
Мне удалось загрузить шейп-файл в python, но я не могу продвинуться дальше.
Это адаптация ответа yosukesabai.
Я хотел убедиться, что искомая точка находится в той же системе проекций, что и шейп-файл, поэтому я добавил для этого код.
Я не мог понять, почему он проводил тест на содержание ply = feat_in.GetGeometryRef()
(в моем тестировании все работало так же хорошо и без него), поэтому я удалил это.
Я также улучшил комментарии, чтобы лучше объяснить, что происходит (насколько я понимаю).
#!/usr/bin/python
import ogr
from IPython import embed
import sys
drv = ogr.GetDriverByName('ESRI Shapefile') #We will load a shape file
ds_in = drv.Open("MN.shp") #Get the contents of the shape file
lyr_in = ds_in.GetLayer(0) #Get the shape file's first layer
#Put the title of the field you are interested in here
idx_reg = lyr_in.GetLayerDefn().GetFieldIndex("P_Loc_Nm")
#If the latitude/longitude we're going to use is not in the projection
#of the shapefile, then we will get erroneous results.
#The following assumes that the latitude longitude is in WGS84
#This is identified by the number "4326", as in "EPSG:4326"
#We will create a transformation between this and the shapefile's
#project, whatever it may be
geo_ref = lyr_in.GetSpatialRef()
point_ref=ogr.osr.SpatialReference()
point_ref.ImportFromEPSG(4326)
ctran=ogr.osr.CoordinateTransformation(point_ref,geo_ref)
def check(lon, lat):
#Transform incoming longitude/latitude to the shapefile's projection
[lon,lat,z]=ctran.TransformPoint(lon,lat)
#Create a point
pt = ogr.Geometry(ogr.wkbPoint)
pt.SetPoint_2D(0, lon, lat)
#Set up a spatial filter such that the only features we see when we
#loop through "lyr_in" are those which overlap the point defined above
lyr_in.SetSpatialFilter(pt)
#Loop through the overlapped features and display the field of interest
for feat_in in lyr_in:
print lon, lat, feat_in.GetFieldAsString(idx_reg)
#Take command-line input and do all this
check(float(sys.argv[1]),float(sys.argv[2]))
#check(-95,47)
Этот сайт, этот сайт и этого сайта помогли проверить прогноз. EPSG:4326
Другой вариант — использовать Shapely (библиотека Python, основанная на GEOS, движке для PostGIS) и Fiona (которая в основном предназначена для чтения/записи файлов):
import fiona
import shapely
with fiona.open("path/to/shapefile.shp") as fiona_collection:
# In this case, we'll assume the shapefile only has one record/layer (e.g., the shapefile
# is just for the borders of a single country, etc.).
shapefile_record = fiona_collection.next()
# Use Shapely to create the polygon
shape = shapely.geometry.asShape( shapefile_record['geometry'] )
point = shapely.geometry.Point(32.398516, -39.754028) # longitude, latitude
# Alternative: if point.within(shape)
if shape.contains(point):
print "Found shape for point."
Обратите внимание, что выполнение тестов «точка в полигоне» может быть дорогостоящим, если полигон большой/сложный (например, шейп-файлы для некоторых стран с чрезвычайно неровными береговыми линиями). В некоторых случаях может помочь использование ограничительных рамок, чтобы быстро исключить вещи, прежде чем выполнять более интенсивный тест:
minx, miny, maxx, maxy = shape.bounds
bounding_box = shapely.geometry.box(minx, miny, maxx, maxy)
if bounding_box.contains(point):
...
Наконец, имейте в виду, что загрузка и анализ больших/неправильных шейп-файлов занимает некоторое время (к сожалению, хранение таких типов полигонов в памяти также часто обходится дорого).
shapely
/fiona
намного проще, чем использование ogr
. Еще одна вещь, которую я нашел полезной: если вы используете несколько слоев/функций в одном шейп-файле, вы можете просто перебирать элементы в коллекции, использовать метод asShape
и возвращать правильный объект, если он содержит вашу точку. Это также можно использовать для сохранения любых атрибутов (таких как имена и коды), которые могут присутствовать в объекте.
- person Husky; 04.01.2021
Вот простое решение, основанное на pyshp и shapely.
Предположим, что ваш шейп-файл содержит только один полигон (но вы можете легко адаптировать его для нескольких полигонов):
import shapefile
from shapely.geometry import shape, Point
# read your shapefile
r = shapefile.Reader("your_shapefile.shp")
# get the shapes
shapes = r.shapes()
# build a shapely polygon from your shape
polygon = shape(shapes[0])
def check(lon, lat):
# build a shapely point from your geopoint
point = Point(lon, lat)
# the contains function does exactly what you want
return polygon.contains(point)
я сделал почти то же самое, что вы делаете вчера, используя gdal ogr с привязкой к python. Это выглядело так.
import ogr
# load the shape file as a layer
drv = ogr.GetDriverByName('ESRI Shapefile')
ds_in = drv.Open("./shp_reg/satreg_etx12_wgs84.shp")
lyr_in = ds_in.GetLayer(0)
# field index for which i want the data extracted
# ("satreg2" was what i was looking for)
idx_reg = lyr_in.GetLayerDefn().GetFieldIndex("satreg2")
def check(lon, lat):
# create point geometry
pt = ogr.Geometry(ogr.wkbPoint)
pt.SetPoint_2D(0, lon, lat)
lyr_in.SetSpatialFilter(pt)
# go over all the polygons in the layer see if one include the point
for feat_in in lyr_in:
# roughly subsets features, instead of go over everything
ply = feat_in.GetGeometryRef()
# test
if ply.Contains(pt):
# TODO do what you need to do here
print(lon, lat, feat_in.GetFieldAsString(idx_reg))
Один из способов сделать это — прочитать файл формы ESRI с помощью библиотеки OGR http://www.gdal.org/ogr, а затем используйте библиотеку геометрии GEOS http://trac.osgeo.org/geos/ для проверки точки в полигоне. Это требует некоторого программирования на C/C++.
Существует также интерфейс Python для GEOS по адресу http://sgillies.net/blog/14/python-geos-module/ (которым я никогда не пользовался). Может быть, это то, что вы хотите?
Другим решением является использование библиотеки http://geotools.org/. То есть на Яве.
Для этого у меня также есть собственное программное обеспечение Java (которое вы можете загрузить с http://www.mapyrus.org плюс jts.jar
с http://www.vividsolutions.com/products.asp ). Вам нужен только текстовый командный файл inside.mapyrus
, содержащий следующие строки, чтобы проверить, лежит ли точка внутри первого многоугольника в файле формы ESRI:
dataset "shapefile", "us_states.shp"
fetch
print contains(GEOMETRY, -120, 46)
И запустить с:
java -cp mapyrus.jar:jts-1.8.jar org.mapyrus.Mapyrus inside.mapyrus
Будет напечатано 1, если точка находится внутри, иначе 0.
Вы также можете получить хорошие ответы, если разместите этот вопрос на https://gis.stackexchange.com/.
Оформить заказ http://geospatialpython.com/2011/01/point-in-polygon.html и http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html
Если вы хотите узнать, какой полигон (из заполненного ими шейп-файла) содержит заданную точку (и у вас тоже есть куча точек), самый быстрый способ — использовать postgis. На самом деле я реализовал версию на основе fiona, используя ответы здесь, но это было мучительно медленно (сначала я использовал многопроцессорность и проверял ограничивающую рамку). 400 минут обработки = 50к баллов. Используя postgis, это заняло менее 10 секунд. Индексы B-дерева эффективны!
shp2pgsql -s 4326 shapes.shp > shapes.sql
Это создаст файл sql с информацией из шейп-файлов, создаст базу данных с поддержкой postgis и запустит этот sql. Создайте индекс gist в столбце geom. Затем, чтобы найти имя многоугольника:
sql="SELECT name FROM shapes WHERE ST_Contains(geom,ST_SetSRID(ST_MakePoint(%s,%s),4326));"
cur.execute(sql,(x,y))