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

2024-10-01 09:16:26 发布

您现在位置:Python中文网/ 问答频道 /正文

我有一个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)

Tags: 文件theinfromlayerfortransformdataset
3条回答

另一个选择是使用Shaviy(基于GeOS的Python库,PASGIS的引擎)和菲奥娜(基本上用于读/写文件):

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

请注意,如果多边形很大/很复杂(例如,一些海岸线极不规则的国家的形状文件),则进行点对多边形测试可能会很昂贵。在某些情况下,在进行更深入的测试之前,可以使用边界框快速排除问题:

^{pr2}$

最后,请记住,加载和解析大的/不规则的shapefile需要一些时间(不幸的是,这些类型的多边形在内存中的保存通常也很昂贵)。在

您可以使用到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...

签出Python Shapefile Library

这应该给你几何和不同的信息。在

相关问题 更多 >