Alternativas al uso de ArcPy

Question

Parece que uso el paquete ArcPy de ESRI para prácticamente todo mi geoprocesamiento en Python. Para el crédito de ESRI, se trata de un increíble conjunto de herramientas que pueden ayudar a lograr mucho. Sin embargo, también me gustaría crear scripts de geoprocesamiento fuera del dominio de ESRI ArcPy. Por ejemplo, si quiero recortar un raster a un polígono, empezaría con el siguiente script de ESRI :

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "C:/sapyexamples/data"

# Set local variables
inRaster = "elevation"
inMaskData = "mask.shp"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute ExtractByMask
outExtractByMask = ExtractByMask(inRaster, inMaskData)

# Save the output 
outExtractByMask.save("C:/sapyexamples/output/extractmask")  

No estoy seguro de cómo podría realizar la misma tarea mediante programación sin ArcPy. Mis preguntas para los programadores serios por ahí: ¿Qué colección de herramientas de Python utiliza para realizar tareas que los usuarios de ESRI realizarían con el paquete del sitio ArcPy? ¿Por dónde empiezo?

Answer 1

GDAL es la herramienta que hay que utilizar. De hecho toda esa llamada es una línea para gdal_rasterize:

gdal_rasterize -l mask -i -burn -9999 mask.shp elevation.tif

si conocieras el valor sin datos del dem

Para un poco de control de Python:

lyr = 'mask'
shp = 'mask.shp'
dem = 'elevation.tif'
ndv = -9999
p = os.Popen('gdal_rasterize -l %s -i -burn %d %s %s' % (lyr,ndv,shp,dem)

donde sus variables podrían establecerse en python

Para un pitón completo:

from osgeo import gdal, ogr
from osgeo.gdalconst import *
shp = ogr.Open('mask.shp')
lyr = shp.GetLayer('mask')
dem = gdal.Open('elevation.tif', GA_Update)
ndv = dem.GetRasterBand(1).GetNoDataValue()
gdal.RasterizeLayer(dem, 1, lyr, None, ndv) # other options, such as transformer func, creation options...
dem = None

Sólo he echado un vistazo a la sintaxis de la API de C, así que mi sintaxis para python probablemente esté un poco equivocada. Ver gdal_alg.h: http://gdal.org/gdal__alg_8h.html

Answer 2

Encontrará otras preguntas similares en este sitio que hacen la misma pregunta básica y tienen muy buenas referencias. La más parecida (y detallada) es:

• ¿Cuáles son las herramientas/módulos/add-ins de Python cruciales en los SIG?

Otros son:

• Ejemplos de scripts en Python para el geoprocesamiento de shapefiles sin usar arcpy
• Biblioteca de Python puro para operaciones geométricas
• ¿Qué herramientas en Python están disponibles para hacer la distancia de un gran círculo + la creación de una línea?
• Módulo de Python para eliminar las características de SHP (sin el SIG de escritorio instalado)

Answer 3

En muchas de mis investigaciones académicas trabajo con datos LiDAR haciendo análisis de superficies para la geomorfología. Rápidamente descubrí que realizar muchas operaciones con arcpy era muy lento, especialmente en conjuntos de datos grandes. Como resultado empecé a utilizar:

• pyshp para manipular archivos shape y actualizar tablas de atributos
• numpy para gestionar rásteres ASCII y realizar análisis basados en el kernel, como los cálculos de curvatura
• scipy para realizar un análisis estadístico de los resultados y llevar a cabo un ajuste de curvas para las superficies
• matplotlib para trazar gráficos y otros resultados gráficos, como mapas básicos para visualizaciones rápidas

También recomendaría el libro, Modelización cuantitativa de los procesos de la superficie terrestre a cualquier persona que quiera aprender más sobre el análisis de superficies rasterizadas. El libro viene con grandes ejemplos de código en C++, que son mucho más eficientes que las herramientas de ArcGIS. Estos algoritmos también pueden ser portados a Python sin necesidad de nada más complejo que numpy, aunque se ejecutan mucho más rápido en C++.

Answer 4

Un buen punto de partida sería el Biblioteca de abstracción de datos geoespaciales . En realidad se compone de dos bibliotecas: GDAL para manipular datos geoespaciales rasterizados y OGR para manipular datos geoespaciales vectoriales, pero la gente suele llamarla simplemente GDAL.

Hay un geoprocesamiento con Python utilizando SIG de código abierto curso en la Universidad del Estado de Utah. Quizá también quieras comprobarlo.

Answer 5

Creo que las respuestas dadas hasta el momento cubren básicamente todos los paquetes que merecen ser mencionados (especialmente GDAL, OGR, pyshp, NumPy)

Pero también existe la Laboratorio de software GIS y Python que alberga un par de módulos interesantes. Estos son:

• Fiona : La API más ordenada de OGR
• Rtree Índice espacial para Python GIS
• Shapely : Paquete de Python para la manipulación y el análisis de características en el plano cartesiano

Personalmente, he empezado a jugar con GDAL/OGR últimamente y los he encontrado muy impresionantes con respecto a la velocidad y la cobertura de las herramientas de análisis.

Aquí algunos ejemplos de cómo utilizar los métodos (tomados de este excelente fuente que es un muy buen punto de partida):

# To select by attribute:
.SetAttributeFilter("soil = 'clay'")

# To select by location, either:
.SetSpatialFilter(<geom>)   

# or
.SetSpatialFilterRect(<minx>, <miny>, <maxx>, <maxy>)

# DataSource objects have a method `ExecuteSQL(<SQL>)`
.ExecuteSQL("SELECT* FROM sites WHERE soil = 'clay' ORDER BY id DESC")

# Plus all the well known tools, like:

# intersect
poly2.Intersect(<geom_1>)

# disjoint?
<geom>.Disjoint(geom_1)

# touches (on the edge?)
<geom>.Touches(geom_1)

# cross each other?
<geom>.Crosses(geom_1)

# within?
<geom>.Within(geom_1)

#contains?
<geom>.Contains(ptB)

# overlaps?
<geom>.Overlaps(geom_1)

## geoprecessing
<geom>.Union(<geom_1>)
<geom>.Intersection(<geom_1>)
<geom>.Difference(<geom_1>)
<geom>.SymmetricDifference(<geom_1>)

# Buffer (returns a new geometry)
<geom>.Buffer(<distance>)

# Are the geometries equal?
<geom1>.Equal(<geom2>)

# Returns the shortest distance between the two geometries
<geom1>.Distance(<geom2>)

# Returns the geometry's extent as a list (minx, maxx, miny, maxy)
<geom>.GetEnvelope()

Lo bueno de estas herramientas es que son muy flexibles a la hora de aplicarlas. Por ejemplo, escribí mi propia clase CreateGeometry() para crear fácilmente archivos vectoriales desde cero. Si te interesa también puedo publicarlo aquí, aunque creo que está fuera del alcance de la pregunta.