Как использовать python для поиска информации на определенной широте/долготе в шейп-файле ESRI?

У меня есть шейп-файл ESRI (отсюда: 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 gis shapefile
person arkottke    schedule 04.01.2011    source источник
comment
Если точка не найдена ни в одном из features, то последние объекты будут считаться совпадением (и будут напечатаны его поля). Я бы объявил отдельную переменную matched_feature и присвоил ей значение непосредственно перед break, а затем использовал ее для следующего цикла вместо переменной feature   -  person Simon Ndunda    schedule 25.09.2017

Ответы (3)


arrow_upward
2
arrow_downward

Вы можете использовать привязки 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...
person albertov    schedule 04.01.2011
comment
Спасибо, я рассмотрю это решение. - person arkottke; 05.01.2011
comment
У меня возникли некоторые трудности с указанием точки, потому что необходимо преобразовать координаты из широты/долготы в систему отсчета, используемую шейп-файлом. Я предполагаю, что это где-то в наборе инструментов gdal, но я не могу найти его... пока. - person arkottke; 05.01.2011
comment
Посмотрите на osgeo.osr.SpatialReference и geom.Transform() (IIRC). Я обновлю пример завтра (afk) - person albertov; 05.01.2011
comment
Да, я нашел помощь в IRC и провел меня через трансформацию. Я обновил OP, чтобы отразить эту информацию. Спасибо за вашу помощь. - person arkottke; 05.01.2011

arrow_upward
3
arrow_downward

Ознакомьтесь с библиотекой шейп-файлов Python.

Это должно дать вам геометрию и другую информацию.

person milkypostman    schedule 04.01.2011
comment
Да, он предоставляет геометрию, но не дает возможности увидеть, находится ли точка в пределах указанной геометрии. - person arkottke; 05.01.2011

arrow_upward
2
arrow_downward

Другой вариант — использовать Shapely (библиотека Python, основанная на GEOS, движке для PostGIS) и 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."

Обратите внимание, что выполнение тестов «точка в полигоне» может быть дорогостоящим, если полигон большой/сложный (например, шейп-файлы для некоторых стран с чрезвычайно неровными береговыми линиями). В некоторых случаях может помочь использование ограничительных рамок, чтобы быстро исключить вещи, прежде чем выполнять более интенсивный тест:

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

if bounding_box.contains(point):
    ...

Наконец, имейте в виду, что загрузка и анализ больших/неправильных шейп-файлов занимает некоторое время (к сожалению, хранение таких типов полигонов в памяти также часто обходится дорого).

person Clint Harris    schedule 11.09.2013
comment
Спасибо, что указали на fiona, похоже, это отличный пакет. - person arkottke; 29.11.2013