26

ジオポイントが特定のシェープファイルの領域内にあるかどうかを確認するにはどうすればよいですか？

Pythonでシェープファイルをロードできましたが、それ以上取得できません。

4

7 に答える 7

32

もう1つのオプションは、Shapely（PostGISのエンジンであるGEOSに基づくPythonライブラリ）と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):
    ...

最後に、大きな/不規則なシェープファイルのロードと解析には時間がかかることに注意してください（残念ながら、これらのタイプのポリゴンは、メモリに保持するのに費用がかかることがよくあります）。

于 2013-09-11T19:11:49.743 に答える
21

これは、陽介さばいの答えを改作したものです。

検索したポイントがシェープファイルと同じプロジェクションシステムにあることを確認したかったので、そのためのコードを追加しました。

なぜ彼がcontainsテストを行っているのか理解できなかったので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

于 2012-11-17T17:49:28.060 に答える
14

これは、 pyshpとshapelyに基づく簡単なソリューションです。

シェープファイルに含まれるポリゴンは1つだけであると仮定します（ただし、複数のポリゴンに簡単に適応できます）。

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

昨日、Pythonバインディングでgdalのogrを使用して、ほぼ正確に実行しました。こんな感じでした。

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

これを行う1つの方法は、OGRライブラリhttp://www.gdal.org/ogrを使用してESRIシェープファイルを読み取り、次にGEOSジオメトリライブラリhttp://trac.osgeo.org/geos/を使用してポイントを実行することです。 -ポリゴン内テスト。これには、C /C++プログラミングが必要です。

http://sgillies.net/blog/14/python-geos-module/にGEOSへのPythonインターフェースもあります（これは私が使用したことはありません）。多分それはあなたが望むものですか？

別の解決策は、 http： //geotools.org/ライブラリを使用することです。それはJavaです。

これを行うための独自のJavaソフトウェアもあります（http://www.mapyrus.orgおよびhttp://www.vividsolutions.com/products.aspjts.jarからダウンロードできます）。ESRIシェープファイルの最初のポリゴンの内側にポイントがあるかどうかを確認するには、次の行を含むテキストコマンドファイルのみが必要です。inside.mapyrus

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 に答える
1

チェックアウトhttp://geospatialpython.com/2011/01/point-in-polygon.htmlおよびhttp://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html

于 2011-10-24T15:47:31.807 に答える
0

（それらでいっぱいのシェープファイルから）どのポリゴンに特定のポイントが含まれているかを知りたい場合（そしてポイントがたくさんある場合）、最速の方法はpostgisを使用することです。ここでの回答を使用して、実際にフィオナベースのバージョンを実装しましたが、非常に低速でした（最初にマルチプロセッシングを使用し、バウンディングボックスをチェックしていました）。400分の処理=50kポイント。postgisを使用すると、10秒もかかりませんでした。Bツリーインデックスは効率的です！

shp2pgsql -s 4326 shapes.shp > shapes.sql

これにより、シェープファイルからの情報を含むSQLファイルが生成され、postgisをサポートするデータベースが作成され、そのSQLが実行されます。geom列に要点インデックスを作成します。次に、ポリゴンの名前を見つけるには：

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 に答える