Проверьте, находится ли геопункт с широтой и долготой в шейп-файле - PullRequest
Новая
       52

Проверьте, находится ли геопункт с широтой и долготой в шейп-файле

20 голосов
/ 22 октября 2011

Как я могу проверить, находится ли геопункт в пределах области данного шейп-файла?

Мне удалось загрузить шейп-файл в python, но я не могу продолжить.

Ответы [ 7 ]

30 голосов
/ 11 сентября 2013

Другой вариант - использовать 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):
    ...

И наконец, имейте в виду, что загрузка и анализ больших / неправильных шейп-файлов занимает некоторое время (к сожалению, такие типы полигонов зачастую слишком дороги для хранения в памяти).

16 голосов
/ 17 ноября 2012

Это адаптация ответа Йосукесабая.

Я хотел убедиться, что искомая точка находилась в той же проекционной системе, что и шейп-файл, поэтому я добавил код для этого.

Я не мог понять, почему он проводил тест содержимого на 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

8 голосов
/ 09 марта 2017

Вот простое решение, основанное на 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)
2 голосов
/ 27 октября 2011

Я сделал почти то же, что вы делаете вчера, используя ogr gdal с привязкой к 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))
1 голос
/ 24 октября 2011
1 голос
/ 22 октября 2011

Один из способов сделать это - прочитать файл ESRI Shape с помощью 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 Shape: 

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/

0 голосов
/ 07 декабря 2016

Если вы хотите выяснить, какой полигон (из полного в них шейп-файла) содержит данную точку (и у вас также есть куча точек), самый быстрый способ - использовать 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))
...