Pythonを使用して、Esri Shapefileの特定の緯度/経度で情報を検索する方法は？

Question

私はEsri ShapeFileを持っています（ここから： http://pubs.usgs.gov/ds/425/）。特定の緯度/経度で、Pythonを使用してShapeファイル（この場合はサーチフィア材料）から情報を検索したいと考えています。

この問題を解決するための最良の方法は何ですか？

ありがとう。

最終的解決：

#!/usr/bin/python

from osgeo import ogr, osr

dataset = ogr.Open('./USGS_DS_425_SHAPES/Surficial_materials.shp')
layer = dataset.GetLayerByIndex(0)
layer.ResetReading()

# Location for New Orleans: 29.98 N, -90.25 E
point = ogr.CreateGeometryFromWkt("POINT(-90.25 29.98)")

# Transform the point into the specified coordinate system from WGS84
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(4326)
coordTransform = osr.CoordinateTransformation(
        spatialRef, layer.GetSpatialRef())

point.Transform(coordTransform)

for feature in layer:
    if feature.GetGeometryRef().Contains(point):
        break

for i in range(feature.GetFieldCount()):
    print feature.GetField(i)

Answer 1

Python Bindingsをに使用できます GDAL/OGR ツールキット。これが例です：

from osgeo import ogr

ds = ogr.Open("somelayer.shp")
lyr = ds.GetLayerByName("somelayer")
lyr.ResetReading()

point = ogr.CreateGeometryFromWkt("POINT(4 5)")

for feat in lyr:
    geom = feat.GetGeometryRef()
    if geom.Contains(point):
        sm = feat.GetField(feat.GetFieldIndex("surface_material"))
        # do stuff...

Answer 2

チェックアウト Python ShapeFileライブラリ

これにより、ジオメトリとさまざまな情報が得られるはずです。

Answer 3

別のオプションは、Shapely（GEOSに基づくPythonライブラリ、ポストGISのエンジン）と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):
    ...

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