Как использовать Python для поиска информации в определенной широте/долготе в форме ESRI?
Вопрос
У меня есть файл esri Shapefile (отсюда: 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)
Решение
Вы можете использовать привязки Python в GDAL/OGR Инструментарий. Вот пример:
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...
Другие советы
Оформить Библиотека питона Shapefile
Это должно дать вам геометрию и различную информацию.
Другой вариант - использовать Shapely (библиотека Python на основе GEO, двигателя для постгиса) и Fiona (которая в основном предназначена для чтения/написания файлов):
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."
Обратите внимание, что проведение испытаний на точке полигона может быть дорого, если многоугольник большой/сложный (например, файлы Shape для некоторых стран с чрезвычайно нерегулярными береговыми линиями). В некоторых случаях он может помочь использовать ограничительные ящики, чтобы быстро исключить ситуацию, прежде чем пройти более интенсивный тест:
minx, miny, maxx, maxy = shape.bounds
bounding_box = shapely.geometry.box(minx, miny, maxx, maxy)
if bounding_box.contains(point):
...
Наконец, имейте в виду, что для загрузки и разбора больших/нерегулярных форм для загрузки и разбора больших/нерегулярных форм (к сожалению, эти типы многоугольников часто бывают дорогими, чтобы держать в памяти).