Come usare python per cercare informazioni a latitudine/longitudine specifiche in uno shape file ESRI?

Question

Ho uno shapefile ESRI (da qui: http://pubs.usgs.gov/ds/425/). Sto cercando di usare python per cercare informazioni dal file shape (materiale informativo in questo caso) ad una data latitudine/longitudine.

Qual è il modo migliore per risolvere questo problema?

Grazie.

soluzione finale:

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

È possibile utilizzare le associazioni Python per il toolkit gdal/ogr. Ecco un esempio:

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

Acquista il Python Shapefile Library

Questo dovrebbe dare geometria e informazioni diverse.

Answer 3

Un'altra opzione è quella di utilizzare Shapely (una libreria Python basata su GEOS, il motore per PostGIS) e Fiona (che è fondamentalmente per la lettura/scrittura di file):

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

Si noti che eseguire test point-in-poligono può essere costoso se il poligono è grande/complicato (ad esempio, shapefile per alcuni paesi con coste estremamente irregolari). In alcuni casi può aiutare a utilizzare rettangoli di selezione per escludere rapidamente le cose prima di fare il test più intenso:

minx, miny, maxx, maxy = shape.bounds 
bounding_box = shapely.geometry.box(minx, miny, maxx, maxy) 

if bounding_box.contains(point): 
    ... 

Infine, tenere presente che ci vuole del tempo per caricare e analizzare grandi shapefile/irregolari (purtroppo, quei tipi di poligoni sono spesso costosi da conservare anche in memoria).