blocks|key|533681|text|为此，一种方法是使用OGR库http://www.gdal.org/ogr读取ESRI形状文件，然后使用GEOS几何库http://trac.osgeo.org/geos/执行多边形中的点测试。这需要一些C/C%2B%2B编程。|type|unstyled|depth|inlineStyleRanges|entityRanges|offset|length|data|533682|在http://sgillies.net/blog/14/python-geos-module/上还有一个到GEOS的python接口(我从未用过)。也许这就是你想要的？|533683|另一种解决方案是使用http://geotools.org/库。这是用Java编写的。|533684|我也有自己的Java软件可以做到这一点(您可以从http://www.mapyrus.org下载，还可以从http://www.vividsolutions.com/products.asp下载jts.jar+)。只需包含以下行的文本命令文件inside.mapyrus即可检查点是否位于ESRI形文件中的第一个多边形内：|style|CODE|533685|dataset+"shapefile",+"us_states.shp"
fetch
print+contains(GEOMETRY,+-120,+46)|code-block|syntax|javascript|533686|并使用以下命令运行：|533687|java+-cp+mapyrus.jar:jts-1.8.jar+org.mapyrus.Mapyrus+inside.mapyrus|533688|如果点在内部，它将打印1，否则将打印0。|533689|如果你把这个问题发布在https://gis.stackexchange.com/上，你可能也会得到一些好的答案|533690|entityMap|0|LINK|mutability|MUTABLE|url|http://www.gdal.org/ogr|1|http://trac.osgeo.org/geos/|2|http://sgillies.net/blog/14/python-geos-module/|3|http://geotools.org/|4|http://www.mapyrus.org/|5|http://www.vividsolutions.com/products.asp|6|https://gis.stackexchange.com/^0|E|N|0|1N|R|1|0|1|1B|2|0|A|K|3|0|2P|7|3D|E|O|M|4|1H|16|5|0|0|0|0|0|B|U|6|0^^$0|@$1|2|3|4|5|6|7|1I|8|@]|9|@$A|1J|B|1K|1|1L]|$A|1M|B|1N|1|1O]]|C|$]]|$1|D|3|E|5|6|7|1P|8|@]|9|@$A|1Q|B|1R|1|1S]]|C|$]]|$1|F|3|G|5|6|7|1T|8|@]|9|@$A|1U|B|1V|1|1W]]|C|$]]|$1|H|3|I|5|6|7|1X|8|@$A|1Y|B|1Z|J|K]|$A|20|B|21|J|K]]|9|@$A|22|B|23|1|24]|$A|25|B|26|1|27]]|C|$]]|$1|L|3|M|5|N|7|28|8|@]|9|@]|C|$O|P]]|$1|Q|3|R|5|6|7|29|8|@]|9|@]|C|$]]|$1|S|3|T|5|N|7|2A|8|@]|9|@]|C|$O|P]]|$1|U|3|V|5|6|7|2B|8|@]|9|@]|C|$]]|$1|W|3|X|5|6|7|2C|8|@]|9|@$A|2D|B|2E|1|2F]]|C|$]]|$1|Y|3|-4|5|6|7|2G|8|@]|9|@]|C|$]]]|Z|$10|$5|11|12|13|C|$14|15]]|16|$5|11|12|13|C|$14|17]]|18|$5|11|12|13|C|$14|19]]|1A|$5|11|12|13|C|$14|1B]]|1C|$5|11|12|13|C|$14|1D]]|1E|$5|11|12|13|C|$14|1F]]|1G|$5|11|12|13|C|$14|1H]]]]

One way to do this is to read the ESRI Shape file using the OGR
library <a href="http://www.gdal.org/ogr" rel="nofollow noreferrer">http://www.gdal.org/ogr</a> and then use the GEOS geometry
library <a href="http://trac.osgeo.org/geos/" rel="nofollow noreferrer">http://trac.osgeo.org/geos/</a> to do the point-in-polygon test.
This requires some C/C++ programming.

There is also a python interface to GEOS at <a href="http://sgillies.net/blog/14/python-geos-module/" rel="nofollow noreferrer">http://sgillies.net/blog/14/python-geos-module/</a> (which I have never used). Maybe that is what you want?

Another solution is to use the <a href="http://geotools.org/" rel="nofollow noreferrer">http://geotools.org/</a> library.
That is in Java.

I also have my own Java software to do this (which you can download
from <a href="http://www.mapyrus.org" rel="nofollow noreferrer">http://www.mapyrus.org</a> plus <code>jts.jar</code> from <a href="http://www.vividsolutions.com/products.asp" rel="nofollow noreferrer">http://www.vividsolutions.com/products.asp</a> ). You need only a text command
file <code>inside.mapyrus</code> containing
the following lines to check if a point lays inside the
first polygon in the ESRI Shape file:

<pre><code>dataset "shapefile", "us_states.shp"
fetch
print contains(GEOMETRY, -120, 46)
</code></pre>

And run with:

<pre><code>java -cp mapyrus.jar:jts-1.8.jar org.mapyrus.Mapyrus inside.mapyrus
</code></pre>

It will print a 1 if the point is inside, 0 otherwise.

You might also get some good answers if you post this question on
<a href="https://gis.stackexchange.com/">https://gis.stackexchange.com/</a>

blocks|key|533698|text|查看http://geospatialpython.com/2011/01/point-in-polygon.html和http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html|type|unstyled|depth|inlineStyleRanges|entityRanges|offset|length|data|533699|entityMap|0|LINK|mutability|MUTABLE|url|http://geospatialpython.com/2011/01/point-in-polygon.html|1|http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html^0|2|1L|0|1O|1V|1|0^^$0|@$1|2|3|4|5|6|7|N|8|@]|9|@$A|O|B|P|1|Q]|$A|R|B|S|1|T]]|C|$]]|$1|D|3|-4|5|6|7|U|8|@]|9|@]|C|$]]]|E|$F|$5|G|H|I|C|$J|K]]|L|$5|G|H|I|C|$J|M]]]]

Checkout <a href="http://geospatialpython.com/2011/01/point-in-polygon.html" rel="nofollow">http://geospatialpython.com/2011/01/point-in-polygon.html</a> and <a href="http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html" rel="nofollow">http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html</a>

blocks|key|750259|text|我使用gdal的ogr和python绑定几乎做了你昨天正在做的事情。它看起来像这样。|type|unstyled|depth|inlineStyleRanges|entityRanges|data|750260|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))|code-block|syntax|javascript|750261|entityMap^0|0|0^^$0|@$1|2|3|4|5|6|7|I|8|@]|9|@]|A|$]]|$1|B|3|C|5|D|7|J|8|@]|9|@]|A|$E|F]]|$1|G|3|-4|5|6|7|K|8|@]|9|@]|A|$]]]|H|$]]

i did almost exactly what you are doing yesterday using gdal's ogr with python binding. It looked like this.

<pre><code>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))
</code></pre>

blocks|key|533808|text|这是yosukesabai的答案的改编。|type|unstyled|depth|inlineStyleRanges|entityRanges|data|533809|我想确保我要搜索的点与shapefile位于相同的投影系统中，所以我为它添加了代码。|533810|我不明白为什么他要在ply+=+feat_in.GetGeometryRef()上做一个包含测试(在我的测试中，没有它似乎工作得也很好)，所以我删除了它。|offset|length|style|CODE|533811|我还改进了注释，以便更好地解释正在发生的事情(据我所知)。|533812|#!/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)|code-block|syntax|javascript|533813|This+site、this+site和this+site对投影检查很有帮助。EPSG:4326|533814|entityMap|0|LINK|mutability|MUTABLE|url|http://hackmap.blogspot.de/2008/03/ogr-python-projection.html|1|http://gissolved.blogspot.de/2009/05/projections-and-transformation-2.html|2|http://geoinformaticstutorial.blogspot.de/2012/10/reprojecting-shapefile-with-gdalogr-and.html|3|http://spatialreference.org/ref/epsg/wgs-84/^0|0|0|A|U|0|0|0|0|9|0|A|9|1|K|9|2|13|9|3|0^^$0|@$1|2|3|4|5|6|7|16|8|@]|9|@]|A|$]]|$1|B|3|C|5|6|7|17|8|@]|9|@]|A|$]]|$1|D|3|E|5|6|7|18|8|@$F|19|G|1A|H|I]]|9|@]|A|$]]|$1|J|3|K|5|6|7|1B|8|@]|9|@]|A|$]]|$1|L|3|M|5|N|7|1C|8|@]|9|@]|A|$O|P]]|$1|Q|3|R|5|6|7|1D|8|@]|9|@$F|1E|G|1F|1|1G]|$F|1H|G|1I|1|1J]|$F|1K|G|1L|1|1M]|$F|1N|G|1O|1|1P]]|A|$]]|$1|S|3|-4|5|6|7|1Q|8|@]|9|@]|A|$]]]|T|$U|$5|V|W|X|A|$Y|Z]]|10|$5|V|W|X|A|$Y|11]]|12|$5|V|W|X|A|$Y|13]]|14|$5|V|W|X|A|$Y|15]]]]

This is an adaptation of yosukesabai's answer.

I wanted to ensure that the point I was searching for was in the same projection system as the shapefile, so I've added code for that.

I couldn't understand why he was doing a contains test on <code>ply = feat_in.GetGeometryRef()</code> (in my testing things seemed to work just as well without it), so I removed that.

I've also improved the commenting to better explain what's going on (as I understand it).

<pre><code>#!/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)
</code></pre>

<a href="http://hackmap.blogspot.de/2008/03/ogr-python-projection.html" rel="noreferrer">This site</a>, <a href="http://gissolved.blogspot.de/2009/05/projections-and-transformation-2.html" rel="noreferrer">this site</a>, and <a href="http://geoinformaticstutorial.blogspot.de/2012/10/reprojecting-shapefile-with-gdalogr-and.html" rel="noreferrer">this site</a> were helpful regarding the projection check. <a href="http://spatialreference.org/ref/epsg/wgs-84/" rel="noreferrer">EPSG:4326</a>

blocks|key|750442|text|另一种选择是使用Shapely+(一个基于GEOS的Python库，PostGIS的引擎)和Fiona+(主要用于读/写文件)：|type|unstyled|depth|inlineStyleRanges|entityRanges|data|750443|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."|code-block|syntax|javascript|750444|请注意，如果多边形很大/复杂(例如，一些海岸线极其不规则的国家的shapefile)，执行多边形中点测试的成本可能会很高。在某些情况下，在进行更密集的测试之前，使用边界框快速排除可能会有所帮助：|750445|minx,+miny,+maxx,+maxy+=+shape.bounds
bounding_box+=+shapely.geometry.box(minx,+miny,+maxx,+maxy)

if+bounding_box.contains(point):
++++...|750446|最后，请记住，加载和解析大型/不规则的shapefile需要一些时间(不幸的是，这些类型的多边形在内存中保存的代价也很高)。|750447|entityMap^0|0|0|0|0|0^^$0|@$1|2|3|4|5|6|7|O|8|@]|9|@]|A|$]]|$1|B|3|C|5|D|7|P|8|@]|9|@]|A|$E|F]]|$1|G|3|H|5|6|7|Q|8|@]|9|@]|A|$]]|$1|I|3|J|5|D|7|R|8|@]|9|@]|A|$E|F]]|$1|K|3|L|5|6|7|S|8|@]|9|@]|A|$]]|$1|M|3|-4|5|6|7|T|8|@]|9|@]|A|$]]]|N|$]]

Another option is to use Shapely (a Python library based on GEOS, the engine for PostGIS) and Fiona (which is basically for reading/writing files):

<pre><code>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."
</code></pre>

Note that doing point-in-polygon tests can be expensive if the polygon is large/complicated (e.g., shapefiles for some countries with extremely irregular coastlines). In some cases it can help to use bounding boxes to quickly rule things out before doing the more intensive test:

<pre><code>minx, miny, maxx, maxy = shape.bounds
bounding_box = shapely.geometry.box(minx, miny, maxx, maxy)

if bounding_box.contains(point):
 ...
</code></pre>

Lastly, keep in mind that it takes some time to load and parse large/irregular shapefiles (unfortunately, those types of polygons are often expensive to keep in memory, too).

blocks|key|534044|text|如果你想找出哪个多边形(从一个充满它们的shapefile中)包含一个给定点(并且你也有一堆点)，最快的方法是使用postgis。我实际上实现了一个基于fiona的版本，使用了这里的答案，但它非常慢(我首先使用多进程和检查边界框)。400分钟的处理=+50k分。使用postgis，这只需要不到10秒。B树索引是高效的！|type|unstyled|depth|inlineStyleRanges|entityRanges|data|534045|shp2pgsql+-s+4326+shapes.shp+>+shapes.sql|code-block|syntax|javascript|534046|这将使用shapefile中的信息生成一个sql文件，创建一个支持postgis的数据库并运行该sql。在geom列上创建一个gist索引。然后，要查找多边形的名称：|534047|sql="SELECT+name+FROM+shapes+WHERE+ST_Contains(geom,ST_SetSRID(ST_MakePoint(%25s,%25s),4326));"
cur.execute(sql,(x,y))|534048|entityMap^0|0|0|0|0^^$0|@$1|2|3|4|5|6|7|M|8|@]|9|@]|A|$]]|$1|B|3|C|5|D|7|N|8|@]|9|@]|A|$E|F]]|$1|G|3|H|5|6|7|O|8|@]|9|@]|A|$]]|$1|I|3|J|5|D|7|P|8|@]|9|@]|A|$E|F]]|$1|K|3|-4|5|6|7|Q|8|@]|9|@]|A|$]]]|L|$]]

If you want to find out which polygon (from a shapefile full of them) contains a given point (and you have a bunch of points as well), the fastest way is using postgis. I actually implemented a fiona based version, using the answers here, but it was painfully slow (I was using multiprocessing and checking bounding box first). 400 minutes of processing = 50k points. Using postgis, that took less than 10seconds. B tree indexes are efficient!

<pre><code>shp2pgsql -s 4326 shapes.shp &gt; shapes.sql
</code></pre>

That will generate a sql file with the information from the shapefiles, create a database with postgis support and run that sql. Create a gist index on the geom column. Then, to find the name of the polygon:

<pre><code>sql="SELECT name FROM shapes WHERE ST_Contains(geom,ST_SetSRID(ST_MakePoint(%s,%s),4326));"
cur.execute(sql,(x,y))
</code></pre>

blocks|key|750549|text|这是一个基于pyshp和shapely的简单解决方案。|type|unstyled|depth|inlineStyleRanges|entityRanges|offset|length|data|750550|让我们假设您的shapefile只包含一个多边形(但您可以很容易地适应多个多边形)：|750551|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)|code-block|syntax|javascript|750552|entityMap|0|LINK|mutability|MUTABLE|url|https://pypi.python.org/pypi/pyshp|1|https://pypi.python.org/pypi/Shapely^0|6|5|0|C|7|1|0|0|0^^$0|@$1|2|3|4|5|6|7|U|8|@]|9|@$A|V|B|W|1|X]|$A|Y|B|Z|1|10]]|C|$]]|$1|D|3|E|5|6|7|11|8|@]|9|@]|C|$]]|$1|F|3|G|5|H|7|12|8|@]|9|@]|C|$I|J]]|$1|K|3|-4|5|6|7|13|8|@]|9|@]|C|$]]]|L|$M|$5|N|O|P|C|$Q|R]]|S|$5|N|O|P|C|$Q|T]]]]

Here is a simple solution based on <a href="https://pypi.python.org/pypi/pyshp" rel="noreferrer">pyshp</a> and <a href="https://pypi.python.org/pypi/Shapely" rel="noreferrer">shapely</a>.

Let's assume that your shapefile only contains one polygon (but you can easily adapt for multiple polygons):

<pre><code>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)
</code></pre>

How can I check if a geopoint is within the area of a given shapefile? 

I managed to load a shapefile in python, but can't get any further.

Check if a geopoint with latitude and longitude is within a shapefile

数据库

Python

Java

如何检查地理位置是否在给定shapefile的区域内？我设法用python加载了一个shapefile，但是不能再继续了。

问检查具有纬度和经度的geopoint是否在shapefile中
EN

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

问检查具有纬度和经度的geopoint是否在shapefile中EN

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

问检查具有纬度和经度的geopoint是否在shapefile中
EN