É possível observar o conteúdo do Shapefile usando Python sem uma licença do ArcMap?

40

Gostaria de saber se é possível olhar para o conteúdo de um shapefile usando Python sem ter uma licença do ArcMap. A situação é que você pode criar shapefiles a partir de muitos aplicativos diferentes, não apenas do software ESRI. Eu gostaria de criar um script Python que verifique a referência espacial, o tipo de recurso, os nomes e definições de atributos e o conteúdo dos campos em um shapefile e os compare a um conjunto de valores aceitáveis. Eu gostaria que esse script funcionasse mesmo se a organização não tiver nenhuma licença ESRI. Para fazer algo assim, você precisa usar o ArcPy ou pode cavar em um shapefile sem usar o ArcPy?

Caitlin
fonte
11
Depende de quanto esforço você deseja colocar nele. Existem várias bibliotecas de código aberto que ajudarão (eu gosto do OGR conforme a resposta de Aarons), mas se você realmente deseja controlar (e está preparado para trabalhar para isso) o Shapefile (originalmente por Esri) é um formato aberto, consulte en.wikipedia.org/wiki/Shapefile
Michael Stimson
11
Os arquivos de formato ESRI recentes (nos últimos dois anos) estão ocultos em seu novo formato de geodatabase. Parece que nada pode quebrá-los, exceto o software ARCxxx. Muitos órgãos públicos estão usando para informações públicas ... vergonha.

Respostas:

34

Eu recomendaria familiarizar-se com a API Python GDAL / OGR para trabalhar com dados vetoriais e rasterizados. A maneira mais fácil de começar a usar o GDAL / OGR é através de uma distribuição python como python (x, y) , Anaconda ou OSGeo4W .

Mais detalhes sobre o uso do GDAL para suas tarefas específicas:

Além disso, eu recomendaria o seguinte tutorial da USU para você começar.


Tomando emprestado os exemplos acima, o script a seguir usa as ferramentas do FOSS para executar as seguintes ações:

  1. Verifique a referência espacial
  2. Obter campos e tipos de shapefile
  3. Verifique se as linhas em um campo definido pelo usuário contêm algum valor

# 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()
Aaron
fonte
Obrigado pelo insight @MikeT. A documentação do GDAL / OGR usa o método 'Destroy ()' em todo o livro de receitas. Que problemas você vê com esse método?
Aaron
11
Há situações em que segfaults podem ocorrer quando você usa Destroy () e foi um erro de design expor esse método nas ligações. Um método melhor é desreferenciar objetos GDAL inFeature = None. O livro de receitas GDAL / OGR não faz parte ou foi escrito pela equipe principal do GDAL / OGR.
Mike T
@ MikeT Eu editei o post para incluir seus comentários - obrigado.
Aaron
31

Existem muitos módulos para ler shapefiles no Python, mais antigos que o ArcPy, consulte o Python Package Index (PyPi): shapefiles . Também existem muitos exemplos no GIS SE (procure por [Python] Fiona , por exemplo)

Todos podem ler a geometria, os campos e as projeções.

Mas outros módulos como o PySAL: a Biblioteca de Análise Espacial do Python , Cartopy (que usa pyshp ) ou o Matplotlib Basemap também podem ler shapefiles, entre outras coisas.

O mais fácil de usar é o Fiona , mas se você conhece apenas o ArcPy, use o pyshp , porque o osgeo e o Fiona exigem a instalação da biblioteca GDAL C / C ++, o GeoPandas precisa do módulo Pandas e o PySAL é muito grande (muitos, muitos outros tratamentos)

Se você deseja apenas ler o conteúdo de um shapefile, não precisa de coisas complexas, basta usar o protocolo de interface geográfica (GeoJSON) também implementado no ArcPy ( ArcPy: AsShape )

Com Fiona (como dicionários 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')])}

Com pyshp (como dicionários 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'}

Com osgeo / ogr (como dicionários 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}

Com o GeoPandas (como dataframe do 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)

* observação sobre geopandas Você deve usar versões mais antigas do Fiona e GDAL com ele ou ele não será instalado. GDAL: 1.11.2 Fiona: 1.6.0 Geopandas: 0.1.0.dev-

Existem muitos tutoriais na Web e até livros ( Desenvolvimento Geoespacial de Python , Aprendendo Análise Geoespacial com Python e Geoprocessamento com Python , no prelo)

De maneira mais geral, se você deseja usar o Python sem o ArcPy, veja o mapeamento temático simples do shapefile usando o Python?

gene
fonte
Observe que a página principal de Fiona diz: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 21/09/16
2
Evidentemente, a questão é sobre shapefiles e não rasters. Eles são outros módulos para arquivos rasterizados.
gene
Ótima resposta! Algo para atualizar em 2017?
Michael