Ejemplos de Script de Python para geoprocesamiento shapefiles sin utilizar arcpy

Question

Me gustaría utilizar un script en Python que no se basa en arcpy hacer cosas como consulta un shapefile por atributos, crear nueva capa de la selección y calcular áreas de un polígono y convertir polígonos a puntos.

¿Alguien tiene algún ejemplo de código de otros módulos de Python o bibliotecas? Soy capaz de hacer esto fácilmente utilizando arcpy pero quería explorar otras opciones.

Answer 1

Que extraño, como si la gente de repente se descubrió el poder de Python (sin ArcPy que es solo un módulo de Python entre otros), ver por ejemplo la cuestión de Visualizar shapefile en Python:

• geoespacial de procesamiento en Python tiene una historia muy larga, mucho mayor que el de Arcpy (o arcgisscripting) -> no es "imitar" las capacidades de ArcPy aquí, como dice Pablo, la mayoría ya estaban allí antes de ArcPy.
• la referencia para los módulos de Python es la de Python Package Index (Pypi), y hay una sección dedicada: Tema :: Científica/Ingeniería :: SIG
• usted puede hacer nada con estos módulos y a menudo es más fácil y más rápido que ArcPy porque es puro Python (no cursores...).
• Bien formada, es uno de estos módulos para el procesamiento geoespacial geometrías -> calcular áreas de polígonos y convertir los polígonos a los puntos..
• si desea procesar las capas vectoriales, hay osgeo/ogr, Fiona o Pyshp (y otros, menos usado) -> consulta de un archivo shapefile por atributos, crear nueva capa a partir de la selección, calcular áreas de polígonos, convertir los polígonos a los puntos de
• para el procesamiento de los rásteres, el estándar es de osgeo/gdal
• para el análisis espacial, hay Pysal
• para 3D, se pueden utilizar otros Científicos módulos como numpy o scipy (algoritmos 3D, las redes, pero también las estadísticas, la geoestadística, 2D o 3D)
• Y no hablo de mapnik, matplotlib/mapa base,Geodjango y ...

Usted puede combinar todos (Pysal con shapely, ...) y se mezclan con las de otros Científicos módulos.

Por lo tanto para la secuencia de Comandos de Python ejemplos, búsqueda de Pyshp Fiona, ogr, gdal o bien proporcionado en el sig.stackexchange o de internet (hay muchos ejemplos, no sólo en inglés).)
Uno de ellos en francés (las secuencias de comandos y las cifras son universales !):
- Python: el Uso de capas de vector y ráster en una perspectiva geológica, sin software GIS
otro en inglés:
- SIG con Python, bien formada, y Fiona
y en español
- Determinación de áreas de polígonos irregulares utilizando las coordenadas de los vértices
en el sig.stackexchange
- El perfil de elevación de 10 km a cada lado de una línea
- La actualización de los Atributos mediante Pyshp
- Cómo crear un 3D shapefile de una trama?
- Secuencia de Comandos de Python para la obtención de la diferencia de elevación entre dos puntos
- etc

El guión presentado por Aaron puede ser escrito más simplemente con Fiona que sólo utiliza diccionarios de Python:

import fiona
with fiona.open('sites.shp', 'r') as input:
    with open('hw1a.txt', 'w') as output:
       for pt in input:
           id = pt['properties']['id']
           cover = pt['properties']['cover']
           x = str(point['geometry']['coordinates'][0])
           y = str(point['geometry']['coordinates'][21])
           output.write(id + ' ' + x + ' ' + y+ ' ' + cover + '\n')

y si uso bien formada, además de:

from shapely.geometry import shape
with fiona.open('sites.shp', 'r') as input:
    with open('hw1a.txt', 'w') as output:
       for pt in input:
           id = pt['properties']['id']
           cover = pt['properties']['cover']
           x = str(shape(pt['geometry']).x)
           y = str(shape(pt['geometry']).y)
           output.write(id + ' ' + x + ' ' + y+ ' ' + cover + '\n')

Hay también dos libros:

Python Geoespacial Desarrollo de Eric Westra.

El aprendizaje de Análisis Geoespacial con Python de Joel Lawhead

Python también es utilizado como un lenguaje de secuencias de comandos en otras aplicaciones de SIG como QGIS (Quantum GIS), GRASS GIS, gvSIG o OpenJump o modeladores 3D como Paraview (y Blender también !). Y puedes utilizar la mayoría de los módulos geoespaciales en todos estos de la aplicación (consulte Visualización de QGIS datos con Blender)

Answer 2

Recomiendo altamente el USU sitio de Geoprocesamiento con Python usando Open Source GIS para empezar. Que utilizan principalmente el GDAL/OGR biblioteca a lo largo de los ejercicios. La instalación de GDAL/OGR puede ser un poco un desafío, así que esta entrada de blog pueden ser útiles para usted: la Instalación de GDAL (y OGR) para Python en Windows. También echa un vistazo a las Alternativas al uso de Arcpy en el SIG.SE.

El siguiente opensource de geoprocesamiento ejemplo de secuencia de comandos (desde el USU sitio) es utilizado para la extracción de datos de atributos y escribir en un archivo de texto:

# import modules
import ogr, os, sys

# set the working directory
os.chdir('f:/data/classes/python/data')

# open the output text file for writing
file = open('hw1a.txt', 'w')

# get the shapefile driver
driver = ogr.GetDriverByName('ESRI Shapefile')

# open the data source
datasource = driver.Open('sites.shp', 0)
if datasource is None:
  print 'Could not open file'
  sys.exit(1)

# get the data layer
layer = datasource.GetLayer()

# loop through the features in the layer
feature = layer.GetNextFeature()
while feature:

  # get the attributes
  id = feature.GetFieldAsString('id')
  cover = feature.GetFieldAsString('cover')

  # get the x,y coordinates for the point
  geom = feature.GetGeometryRef()
  x = str(geom.GetX())
  y = str(geom.GetY())

  # write info out to the text file
  file.write(id + ' ' + x + ' ' + y + ' ' + cover + '\n')

  # destroy the feature and get a new one
  feature.Destroy()
  feature = layer.GetNextFeature()

# close the data source and text file
datasource.Destroy()
file.close()

Answer 3

Usted podría estar interesado en GDAL/OGR.

GDAL se utiliza para el procesamiento de los rásteres mientras OGR se utiliza para los vectores. Ambos son bibliotecas de código abierto.

Si usted está buscando para quitar alguna dependencia en ArcPy, puede imitar algunas capacidades mediante la lectura de la información a una matriz y el funcionamiento de su propio cálculos en Python puro.

Yo recientemente lo hizo con la selección de los puntos de un polígono, como se ve aquí. Se utiliza el ray casting algoritmo para determinar si un punto está dentro de un polígono, dadas las coordenadas de los vértices del polígono.

Answer 4

Nunca he usado esto personalmente, pero otros en la oficina, como usar bien: https://pypi.python.org/pypi/Shapely