¿Es posible mirar el contenido de Shapefile usando Python sin una licencia de ArcMap?

40

Me preguntaba si es posible mirar el contenido de un archivo de forma usando Python sin tener una licencia de ArcMap. La situación es que puede crear archivos de forma desde muchas aplicaciones diferentes, no solo desde el software ESRI. Me gustaría crear un script de Python que verifique la referencia espacial, el tipo de entidad, los nombres de atributos y las definiciones, y el contenido de los campos en un archivo shape y los compare con un conjunto de valores aceptables. Me gustaría que este script funcione incluso si la organización no tiene ninguna licencia de ESRI. Para hacer algo como esto, ¿tiene que usar ArcPy o puede excavar en un archivo de forma sin usar ArcPy?

Caitlin
fuente
1
Depende de cuánto esfuerzo quieras poner en él ... hay varias bibliotecas de código abierto que te ayudarán (me gusta OGR según la respuesta de Aarons) pero si realmente quieres controlar (y estás preparado para trabajar para ello) el Shapefile (originalmente por Esri) es un formato abierto, ver en.wikipedia.org/wiki/Shapefile
Michael Stimson
1
Los archivos de forma ESRI recientes (últimos años) están ocultos en su nuevo formato de geodatabase. Parece que nada puede romperlos, excepto el software ARCxxx. Muchas agencias públicas lo están utilizando para información pública ... vergüenza.

Respuestas:

34

Recomendaría familiarizarse con la API Python GDAL / OGR para trabajar con datos vectoriales y ráster. La forma más fácil de comenzar a usar GDAL / OGR es a través de una distribución de python como python (x, y) , Anaconda u OSGeo4W .

Más detalles sobre el uso de GDAL para sus tareas específicas:

Además, recomendaría el siguiente tutorial de USU para comenzar.


Tomando prestado de los ejemplos anteriores, el siguiente script usa herramientas FOSS para realizar las siguientes acciones:

  1. Verifique la referencia espacial
  2. Obtener campos y tipos de archivos de forma
  3. Compruebe si las filas en un campo definido por el usuario contienen algún 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
fuente
Gracias por la información @MikeT. La documentación de GDAL / OGR utiliza el método 'Destroy ()' en todo su libro de cocina. ¿Qué problemas ves con ese método?
Aaron
1
Hay situaciones en las que pueden ocurrir segfaults cuando usas Destroy (), y fue un error de diseño exponer este método en los enlaces. Un mejor método es desreferenciar objetos GDAL como inFeature = None. El libro de cocina de GDAL / OGR no forma parte ni está escrito por el equipo central de GDAL / OGR.
Mike T
@MikeT He editado la publicación para incluir sus comentarios, gracias.
Aaron
31

Hay muchos módulos para leer shapefiles en Python, anteriores a ArcPy, mira el Python Package Index (PyPi): shapefiles . También hay muchos ejemplos en GIS SE (busque [Python] Fiona , por ejemplo)

Todos pueden leer la geometría, los campos y las proyecciones.

Pero otros módulos como PySAL: la Biblioteca de análisis espacial de Python , Cartopy (que usa pyshp ) o Matplotlib Basemap también pueden leer archivos de forma, entre otras cosas.

El más fácil de usar es Fiona , pero si solo conoce ArcPy, use pyshp , porque osgeo y Fiona requieren que se instale la biblioteca GDAL C / C ++, GeoPandas necesita el módulo Pandas y PySAL es demasiado grande (muchos, muchos otros tratamientos)

Si solo desea leer el contenido de un shapefile, no necesita cosas complejas, simplemente use el protocolo de interfaz geográfica (GeoJSON) también implementado en ArcPy ( ArcPy: AsShape )

Con Fiona (como diccionarios de 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')])}

Con pyshp (como diccionarios de 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'}

Con osgeo / ogr (como diccionarios de 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}

Con GeoPandas (como marco de datos de 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)

* nota sobre geopandas Debe usar versiones anteriores de Fiona y GDAL o no se instalará. GDAL: 1.11.2 Fiona: 1.6.0 Geopandas: 0.1.0.dev-

Hay muchos tutoriales en la Web e incluso libros ( Desarrollo Geoespacial de Python , Aprendizaje del Análisis Geoespacial con Python y Geoprocesamiento con Python , en prensa)

En términos más generales, si desea usar Python sin ArcPy, consulte Mapeo temático simple de shapefile usando Python.

gene
fuente
Tenga en cuenta que la página principal de Fiona diceThe 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
Evidentemente, la pregunta es sobre los archivos de forma y no los rásteres. Son otros módulos para archivos ráster.
gen
¡Gran respuesta! ¿Algo para actualizar en 2017?
Michael