26

如何检查地理点是否在给定 shapefile 的区域内？

我设法在python中加载了一个shapefile，但不能再进一步了。

4

7 回答 7

32

另一种选择是使用 Shapely（基于 GEOS 的 Python 库，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."

请注意，如果多边形很大/很复杂（例如，某些海岸线极不规则的国家/地区的 shapefile），则进行多边形中的点测试可能会很昂贵。在某些情况下，在进行更密集的测试之前，它可以帮助使用边界框快速排除事物：

minx, miny, maxx, maxy = shape.bounds
bounding_box = shapely.geometry.box(minx, miny, maxx, maxy)

if bounding_box.contains(point):
    ...

最后，请记住，加载和解析大型/不规则形状文件需要一些时间（不幸的是，这些类型的多边形在内存中通常也很昂贵）。

于 2013-09-11T19:11:49.743 回答
21

这是对 yosukesabai 答案的改编。

我想确保我正在搜索的点与 shapefile 在同一个投影系统中，所以我为此添加了代码。

我不明白他为什么要进行包含测试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)

This site , this site , and this site对投影检查很有帮助。每股收益：4326

于 2012-11-17T17:49:28.060 回答
14

这是一个基于pyshp和shapely的简单解决方案。

假设您的 shapefile 仅包含一个多边形（但您可以轻松适应多个多边形）：

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)
于 2017-03-09T18:08:55.370 回答
2

我使用 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))
于 2011-10-26T20:33:27.330 回答
1

一种方法是使用 OGR 库http://www.gdal.org/ogr读取 ESRI Shape 文件，然后使用 GEOS 几何库http://trac.osgeo.org/geos/来完成这一点-in-polygon 测试。这需要一些 C/C++ 编程。

在http://sgillies.net/blog/14/python-geos-module/也有一个到 GEOS 的 python 接口（我从未使用过）。也许这就是你想要的？

另一种解决方案是使用http://geotools.org/库。那是在Java中。

我也有自己的 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/上发布此问题，您也可能会得到一些好的答案

于 2011-10-22T17:40:16.560 回答
0

如果你想找出哪个多边形（从一个充满它们的 shapefile）包含一个给定的点（你也有一堆点），最快的方法是使用 postgis。我实际上实现了一个基于 fiona 的版本，使用这里的答案，但它非常慢（我首先使用多处理和检查边界框）。400 分钟的处理时间 = 50k 点。使用 postgis，不到 10 秒。B树索引是有效的！

shp2pgsql -s 4326 shapes.shp > shapes.sql

这将使用 shapefile 中的信息生成一个 sql 文件，创建一个支持 postgis 的数据库并运行该 sql。在 geom 列上创建一个 gist 索引。然后，要找到多边形的名称：

sql="SELECT name FROM shapes WHERE ST_Contains(geom,ST_SetSRID(ST_MakePoint(%s,%s),4326));"
cur.execute(sql,(x,y))
于 2016-12-07T17:25:09.320 回答