¿Cómo usar Python para buscar información en latitud/longitud específica en un archivo shape de ESRI?

Question

Tengo un shapefile de ESRI (desde aquí: http://pubs.usgs.gov/ds/425/). Estoy buscando usar Python para buscar información del archivo de formas (material superficial en este caso) en una latitud/longitud determinada.

¿Cuál es la mejor manera de resolver este problema?

Gracias.

solución final:

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

Puede usar los enlaces de python al kit de herramientas gdal/ogr. He aquí un ejemplo:

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

Pedido del Python Shapefile Library

Esto debe darle la geometría y la información diferente.

Answer 3

Otra opción es utilizar bien proporcionado (una biblioteca de Python basado en GEOS, el motor de PostGIS) y Fiona (que es básicamente para la lectura/escritura de archivos):

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

Tenga en cuenta que hacer pruebas de punto en un polígono puede ser costoso si el polígono es grande/complicado (por ejemplo, shapefiles para algunos países con costas extremadamente irregulares). En algunos casos puede ayudar el uso de cuadros delimitadores para descartar rápidamente las cosas antes de hacer la prueba más intensivo:

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

if bounding_box.contains(point): 
    ... 

Por último, tenga en cuenta que se necesita un cierto tiempo para cargar y analizar archivos de formas grandes/irregulares (por desgracia, esos tipos de polígonos a menudo son caros de guardar en la memoria, también).