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

39

Me pregunté si es posible mirar el contenido de un shapefile usando Python sin tener una licencia de ArcMap. La situación es que puede crear shapefiles desde muchas aplicaciones diferentes, no solo desde el software ESRI. Me gustaría crear una secuencia de comandos de Python que compruebe la referencia espacial, el tipo de entidad, los nombres de los atributos y las definiciones, y el contenido de los campos en un shapefile y los compare con un conjunto de valores aceptables. Me gustaría que esta secuencia de comandos funcione incluso si la organización no tiene licencias ESRI. Para hacer algo como esto, ¿tienes que usar ArcPy o puedes profundizar en un shapefile sin usar ArcPy?

    
pregunta Caitlin 04.05.2015 - 19:30

3 respuestas

32

Recomendaría familiarizarse con la API de 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 , o 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 prestados los ejemplos anteriores, la siguiente secuencia de comandos usa herramientas de software libre para realizar las siguientes acciones:

  1. Comprueba la referencia espacial
  2. Obtener campos y tipos de shapefile
  3. Comprueba 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()
    
respondido por el Aaron 04.05.2015 - 19:36
29

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

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

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

Si solo quieres leer el contenido de un shapefile, no necesitas cosas complejas, simplemente usa la interfaz geo protocolo (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 Tienes que usar versiones anteriores 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 ( Python Geospatial Development , Análisis geoespacial de aprendizaje con Python y Geoprocesamiento con Python , en prensa)

Más generalmente, si desea usar Python sin ArcPy, consulte Mapeo temático simple de shapefile usando Python?

    
respondido por el gene 04.05.2015 - 21:19
10

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

La biblioteca de archivos de formas de Python (pyshp)

GeoPandas

Si está interesado en otras bibliotecas, esta publicación en bibliotecas geoespaciales de Python esenciales Es un buen lugar para mirar.

    
respondido por el juturna 04.05.2015 - 19:34

Lea otras preguntas en las etiquetas