Можно ли посмотреть содержимое Shapefile, используя Python без лицензии ArcMap?

40

Я задавался вопросом, можно ли посмотреть содержимое шейп-файла с использованием Python без лицензии ArcMap. Ситуация такова, что вы можете создавать шейп-файлы из множества различных приложений, а не только из программного обеспечения ESRI. Я хотел бы создать скрипт Python, который проверяет пространственную привязку, тип объекта, имена и определения атрибутов, а также содержимое полей в шейп-файле и сравнивает их с набором допустимых значений. Мне бы хотелось, чтобы этот сценарий работал, даже если у организации нет лицензий ESRI. Чтобы сделать что-то вроде этого, вам нужно использовать ArcPy или вы можете копать шейп-файл без использования ArcPy?

Кейтлин
источник
1
Это зависит от того, сколько усилий вы хотите приложить к этому ... есть несколько библиотек с открытым исходным кодом, которые помогут (мне нравится OGR согласно ответу Aarons), но если вы действительно хотите контролировать (и готовы работать для него) Shapefile (первоначально Esri) - открытый формат, см. en.wikipedia.org/wiki/Shapefile
Майкл Стимсон,
1
Последние (последние пару лет) шейп-файлы ESRI скрыты в их новом формате базы геоданных. Кажется, что ничто не может сломать их, кроме программного обеспечения ARCxxx. Многие государственные органы используют его для публичной информации ... позор.

Ответы:

34

Я бы порекомендовал ознакомиться с Python GDAL / OGR API для работы как с векторными, так и с растровыми данными. Самый простой способ начать использовать GDAL / OGR - через дистрибутив python, такой как python (x, y) , Anaconda или OSGeo4W .

Более подробная информация об использовании GDAL для ваших конкретных задач:

Кроме того, я бы порекомендовал вам следующее руководство от USU.


Исходя из приведенных выше примеров, следующий скрипт использует инструменты FOSS для выполнения следующих действий:

  1. Проверьте пространственную привязку
  2. Получить поля и типы шейп-файлов
  3. Проверьте, содержат ли строки в определенном пользователем поле какое-либо значение

# Import the necessary modules
from  osgeo import ogr, osr

driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.Open(r'C:\your\shapefile.shp')

# Get Projection from layer
layer = shp.GetLayer()
spatialRef = layer.GetSpatialRef()
print spatialRef

# Get Shapefile Fields and Types
layerDefinition = layer.GetLayerDefn()

print "Name  -  Type  Width  Precision"
for i in range(layerDefinition.GetFieldCount()):
    fieldName =  layerDefinition.GetFieldDefn(i).GetName()
    fieldTypeCode = layerDefinition.GetFieldDefn(i).GetType()
    fieldType = layerDefinition.GetFieldDefn(i).GetFieldTypeName(fieldTypeCode)
    fieldWidth = layerDefinition.GetFieldDefn(i).GetWidth()
    GetPrecision = layerDefinition.GetFieldDefn(i).GetPrecision()
    print fieldName + " - " + fieldType+ " " + str(fieldWidth) + " " + str(GetPrecision)

# Check if rows in attribute table meet some condition
inFeature = layer.GetNextFeature()
while inFeature:

    # get the cover attribute for the input feature
    cover = inFeature.GetField('cover')

    # check to see if cover == grass
    if cover == 'trees':
        print "Do some action..."

    # destroy the input feature and get a new one
    inFeature = None
    inFeature = inLayer.GetNextFeature()
Аарон
источник
Спасибо за понимание @MikeT. Документация GDAL / OGR использует метод Destroy () во всей своей кулинарной книге. Какие проблемы вы видите с этим методом?
Аарон
1
Существуют ситуации, когда при использовании Destroy () могут происходить ошибки segfa, и было ошибкой проектирования выставлять этот метод в привязках. Лучший способ - разыменовать объекты GDAL, как inFeature = None. Кулинарная книга GDAL / OGR не является частью или написана основной командой GDAL / OGR.
Майк Т
@MikeT Я отредактировал пост, чтобы включить ваши комментарии - спасибо.
Аарон
31

Существует много модулей для чтения шейп-файлов в Python, более старых, чем ArcPy, посмотрите индекс пакетов Python (PyPi): шейп-файлы . Есть также много примеров в GIS SE (например, поиск [Python] Fiona )

Все могут прочитать геометрию, поля и проекции.

  • Чем старше OSGeo (GDAL / OGR) , взгляд на Python GDAL / OGR Cookbook , например ,
  • другое решение - Fiona , также основанное на GDAL / OGR, но более простое в использовании (со словарями Python: формат GeoJSON).
  • pyshp (shapefile) - это чистое решение Python
  • GeoPandas используют Fiona для чтения / записи шейп-файлов и Pandas для инструментов анализа данных. Вы должны знать Панд, чтобы использовать это.

Но другие модули, такие как PySAL: библиотека пространственного анализа Python , Cartopy (которая использует pyshp ) или Matplotlib Basemap , помимо прочего, также могут читать шейп-файлы.

Самым простым в использовании является Fiona , но если вы знаете только ArcPy, используйте pyshp , потому что osgeo и Fiona требуют, чтобы была установлена ​​библиотека GDAL C / C ++, для GeoPandas требуется модуль Pandas, а PySAL слишком большой (многие, многие другие процедуры)

Если вы хотите , чтобы прочитать содержание шейпа, вам не нужны сложные вещи, просто использовать гео интерфейс протокол (GeoJSON) также реализован в ArcPy ( ArcPy: AsShape )

С Фионой (как словари Python):

import fiona
with fiona.open('a_shape.shp') as shp:
     # schema of the shapefile
     print shp.schema
     {'geometry': 'Point', 'properties': OrderedDict([(u'DIP', 'int:2'), (u'DIP_DIR', 'int:3'), (u'TYPE', 'str:10')])}
     # projection
     print shp.crs
     {u'lon_0': 4.367486666666666, u'ellps': u'intl', u'y_0': 5400088.438, u'no_defs': True, u'proj': u'lcc', u'x_0': 150000.013, u'units': u'm', u'lat_2': 49.8333339, u'lat_1': 51.16666723333333, u'lat_0': 90}
     for feature in shp:
        print feature              
{'geometry': {'type': 'Point', 'coordinates': (272070.600041, 155389.38792)}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'DIP', 30), (u'DIP_DIR', 130), (u'TYPE', u'incl')])}
{'geometry': {'type': 'Point', 'coordinates': (271066.032148, 154475.631377)}, 'type': 'Feature', 'id': '1', 'properties': OrderedDict([(u'DIP', 55), (u'DIP_DIR', 145), (u'TYPE', u'incl')])}
{'geometry': {'type': 'Point', 'coordinates': (273481.498868, 153923.492988)}, 'type': 'Feature', 'id': '2', 'properties': OrderedDict([(u'DIP', 40), (u'DIP_DIR', 155), (u'TYPE', u'incl')])}

С pyshp (как словари Python)

import shapefile
reader= shapefile.Reader("a_shape.shp")
# schema of the shapefile
print dict((d[0],d[1:]) for d in reader.fields[1:])
{'DIP_DIR': ['N', 3, 0], 'DIP': ['N', 2, 0], 'TYPE': ['C', 10, 0]}
fields = [field[0] for field in reader.fields[1:]]
for feature in reader.shapeRecords():
    geom = feature.shape.__geo_interface__
    atr = dict(zip(fields, feature.record))
    print geom, atr
{'type': 'Point', 'coordinates': (272070.600041, 155389.38792)} {'DIP_DIR': 130, 'DIP': 30, 'TYPE': 'incl'}
{'type': 'Point', 'coordinates': (271066.032148, 154475.631377)} {'DIP_DIR': 145, 'DIP': 55, 'TYPE': 'incl'}
{'type': 'Point', 'coordinates': (273481.498868, 153923.492988)} {'DIP_DIR': 155, 'DIP': 40, 'TYPE': 'incl'}

С osgeo / ogr (как словари Python)

from osgeo import ogr
reader = ogr.Open("a_shape.shp")
layer = reader.GetLayer(0)
for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    print feature.ExportToJson()
{"geometry": {"type": "Point", "coordinates": [272070.60004, 155389.38792]}, "type": "Feature", "properties": {"DIP_DIR": 130, "DIP": 30, "TYPE": "incl"}, "id": 0}
{"geometry": {"type": "Point", "coordinates": [271066.032148, 154475.631377]}, "type": "Feature", "properties": {"DIP_DIR": 145, "DIP": 55, "TYPE": "incl"}, "id": 1}
{"geometry": {"type": "Point", "coordinates": [273481.49887, 153923.492988]}, "type": "Feature", "properties": {"DIP_DIR": 155, "DIP": 40, "TYPE": "incl"}, "id": 2}

С GeoPandas (в качестве кадра данных Pandas)

import geopandas as gp
shp = gp.GeoDataFrame.from_file('a_shape.shp')
print shp
        DIP_DIR    DIP  TYPE                       geometry
0         130       30  incl          POINT (272070.600041 155389.38792)
1         145       55  incl          POINT (271066.032148 154475.631377)
2         155       40  incl          POINT (273481.498868 153923.492988)

* примечание к геопандам. Вы должны использовать более старые версии Fiona и GDAL с ним, иначе он не будет установлен. GDAL: 1.11.2 Фиона: 1.6.0 Геопанды: 0.1.0.dev-

В Интернете есть много учебных пособий и даже книг ( разработка геопространственного анализа Python , изучение геопространственного анализа с использованием Python и геообработка с использованием Python , в печати)

В более общем плане, если вы хотите использовать Python без ArcPy, посмотрите на Простое тематическое отображение шейп-файла с использованием Python?

ген
источник
Обратите внимание, что на главной странице Фионы написаноThe kinds of data in GIS are roughly divided into rasters representing continuous scalar fields (land surface temperature or elevation, for example) and vectors representing discrete entities like roads and administrative boundaries. Fiona is concerned exclusively with the latter
Mawg
2
Очевидно, вопрос касается шейп-файлов, а не растров. Это другие модули для растровых файлов.
ген
Отличный ответ! Что-нибудь обновить в 2017 году?
Майкл
11

Помимо ArcPy, есть геопространственные библиотеки Python, которые дадут вам эти возможности. Вот два примера:

Библиотека шейп-файлов Python (pyshp)

GeoPandas

Если вы интересуетесь другими библиотеками, этот пост о важных библиотеках Python Geospatial будет хорошим местом для поиска.

Ютурны
источник
1
+1 для пышп, супер прост и понятен в использовании. Похоже, все, что хочет ОП, будет учтено.
mr.adam