Jak używać Pythona do wyszukiwania informacji o określonej szerokości/długości geograficznej w pliku kształtu ESRI?

Question

Mam plik kształtu ESRI (tutaj: http://pubs.usgs.gov/ds/425/). Chciałbym użyć Pythona do wyszukiwania informacji z pliku kształtu (materiał powierzchniowy w tym przypadku) na danej szerokości/długości geograficznej.

Jaki jest najlepszy sposób rozwiązania tego problemu?

Dzięki.

rozwiązanie końcowa:

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

Możesz użyć powiązań Pythona z zestawem narzędzi gdal/ogr. Oto przykład:

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

kasie Python Shapefile Library

To powinno dać geometrii i różnych informacji.

Answer 3

Inną opcją jest użycie foremny (biblioteka Pythona w oparciu o GEOS, silnik dla PostGIS) i Fiona (co jest w zasadzie dla odczytu/zapisu plików):

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

Zauważ, że wykonywanie testów punkt-w-poligonach może być kosztowne, jeśli wielokąt jest duży/skomplikowany (np. Shapefile dla niektórych krajów o wyjątkowo nieregularnych liniach brzegowych). W niektórych przypadkach może pomóc użycie obwiednie szybko wykluczyć rzeczy przed wykonaniem bardziej intensywne Test:

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

if bounding_box.contains(point): 
    ... 

Wreszcie, należy pamiętać, że potrzeba trochę czasu, aby załadować i analizować duże/nieregularne Shapefiles (niestety, te typy wielokątów są często kosztowne w utrzymaniu w pamięci).