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

Question

Me preguntaba si es posible mirar el contenido de un shapefile usando Python sin tener una licencia de ArcMap. La situación es que se pueden crear shapefiles desde muchas aplicaciones diferentes, no sólo desde el software de ESRI. Me gustaría crear un script en Python que compruebe la referencia espacial, el tipo de característica, los nombres y definiciones de los atributos y el contenido de los campos de un shapefile y los compare con un conjunto de valores aceptables. Me gustaría que este script funcionara incluso si la organización no tiene ninguna licencia de ESRI. Para hacer algo así, ¿hay que usar ArcPy o se puede indagar en un shapefile sin usar ArcPy?

Answer 1

Recomiendo que se familiarice con el programa Python GDAL/OGR API para trabajar con datos vectoriales y rasterizados. La forma más fácil de empezar a utilizar GDAL/OGR es a través de una distribución de python como python(x,y) , Anaconda o OSGeo4W .

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

• Obtener campos y tipos de Shapefile
• Obtenga la proyección

Además, le recomiendo el siguiente tutorial de la USU para empezar.

• Geoprocesamiento con Python utilizando SIG de código abierto

---

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

1. Comprobar la referencia espacial
2. Obtener campos y tipos de shapefile
3. Comprobar si las filas de 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()

Answer 2

Hay muchos módulos para leer shapefiles en Python, más antiguos que ArcPy, mira el Índice de paquetes de Python (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.

• El mayor es osgeo (GDAL/OGR) , mira el Libro de cocina GDAL/OGR de Python por ejemplo
• otra solución es Fiona También basado en GDAL/OGR, pero más fácil de usar (con diccionarios de Python: formato GeoJSON).
• pyshp (shapefile) es una solución puramente Python
• GeoPandas utilice Fiona para leer/escribir los shapefiles y Pandas para las herramientas de análisis de datos. Es necesario saber Pandas para usarlo.

Pero otros módulos como PySAL: la biblioteca de análisis espacial de Python , Cartopy (que utilizan pyshp ) o Mapa base Matplotlib también puede leer shapefiles, entre otras cosas.

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

Si sólo quieres leer el contenido de un shapefile, no necesitas cosas complejas, simplemente utiliza la función interfaz geográfica (GeoJSON) también implementado en ArcPy ( ArcPy: AsShape )

Con Fiona (como los 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 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 Tienes que usar versiones antiguas de Fiona y GDAL con él 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 en Python , Aprendizaje del análisis geoespacial con Python y Geoprocesamiento con Python , en prensa)

De manera más general, si quiere utilizar Python sin ArcPy, mire ¿Mapeo temático simple de shapefile usando Python?

Answer 3

Hay bibliotecas geoespaciales de Python, además de ArcPy, que le darán estas habilidades. Aquí hay dos ejemplos:

La biblioteca Python Shapefile (pyshp)

GeoPandas

Si está interesado en otras bibliotecas, este post sobre las bibliotecas geoespaciales esenciales de Python es un buen lugar para buscar.