如何使用python在ESRI形状文件中查找特定纬度/经度的信息？

Question

我有一个ESRI形状文件（从这里：http://pubs.usgs.gov/ds/425/）。我希望使用python从给定纬度/经度的形状文件（本例中为表面材质）中查找信息。在

解决这个问题最好的办法是什么？在

谢谢。在

最终解决方案：

#!/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

您可以使用到gdal/ogr工具箱的python绑定。下面是一个例子：

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...