Wie verwendet man Python, um Informationen in einem ESRI-Shapefile nach bestimmten Längen- und Breitengraden zu suchen?

Question

Ich habe ein ESRI Shapefile (von hier: http://pubs.usgs.gov/ds/425/). Ich suche nach Python, um Informationen aus der Shape-Datei (Surficial-Material in diesem Fall) bei einer gegebenen Länge/Breite zu suchen.

Was ist der beste Weg, um dieses Problem zu lösen?

Danke.

Endlösung:

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

Sie können die Python-Bindungen zum gdal/ogr Toolkit verwenden. Hier ein Beispiel:

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

Kasse die Python Shapefile Library

Dies sollte Ihnen und verschiedene Rahmen Geometrie geben.

Answer 3

Eine weitere Option ist Shapely zu verwenden (eine Bibliothek Python basierend auf GEOS, der Motor für PostGIS) und Fiona (das ist im Grunde zum Lesen/Schreiben von Dateien):

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

Beachten Sie, dass Punkt-in-Polygon-Tests teuer sein können, wenn das Polygon groß oder kompliziert ist (z. B. Shapefiles für einige Länder mit extrem unregelmäßigen Küstenlinien). In einigen Fällen kann es helfen, Begrenzungsrahmen zu verwenden, um schnell die Dinge ausschließen, bevor Sie den intensiveren Test machen:

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

if bounding_box.contains(point): 
    ... 

Schließlich bedenken Sie, dass es einige Zeit, groß/unregelmäßige Shape-Dateien zu laden und zu analysieren, nimmt (leider, diese Arten von Polygonen sind oft auch teuer im Gedächtnis zu behalten).