Gabriel Gaona
mié, 08 de enero, 2020
El uso de Python en sistemas de información geográfica requiere por lo menos las siguientes bases de conocimiento:
Aqui dejare que los expertos (La gente de Python Ecuador) les expliquen los detalles de la sintaxis de este lenguaje.
Sistema de
Información
Geográfica
Componentes de un SIG
El problema fundamental que intenta resolver una herramienta de SIG es el analisis de la información con relacion al espacio.
El espacio es la superficie de la tierra en donde se produce la información.
La información puede referirse a objetos o a procesos que se desarrollen en esa superficie.
Modelos de datos en SIG
La información es organizada en estructuras especiales conocidas como capas ~y estas pueden ser almacenadas en ficheros o bases de datos~
Dentro de una capa se encuentran las entidades. Cada entidad posee atributos intrínsecos y agregados
Toda la información es gestionada a diferentes niveles.
Importante! diferenciar entre la gestión información a nivel de software y la gestión a nivel de datos
El uso principal de python en Sistemas de información geográfica está relacionado con “Scripting”.
Recientemente se han introducido intérpretes python para la ampliación de algoritmos y herramientas desde la interfaz de los SIG. Algunos permiten combinar Python con otros lenguajes (ej. bash, C++, Fortan, java, …)
Software_GIS | API | Link |
---|---|---|
ArcGIS 10.x | ArcPy | http://desktop.arcgis.com/en/arcmap/10.3/analyze/arcpy |
QGis 2.x | PyQgis | http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/ |
gvSIG 2.x | gvsig.py | http://docs.gvsig.org/plone/projects/gvsig-desktop/docs/user/gvsig-desktop-2-0-scripting/ |
GRASS-GIS | PyGrass | https://grass.osgeo.org/grass70/manuals/libpython/pygrass_index.html |
SAGA-GIS | SAGA-Pyhton | http://saga-python.readthedocs.org |
Lib. indep. | Gdal/OGR | http://www.gdal.org/ |
Lib. indep. | PySAL | https://pysal.readthedocs.org/en/latest/ |
El problema a resolver es la visualización de estadísticas a partir de los atributos de las entidades que se muestran en la vista.
El ejemplo fue tomado del taller de Scripting con gvSIG de las 11as. Jornadas Internacionales ()
Necesitamos el módulo de scripting de gvSIG que se importará con las siguientes líneas:
from gvsig import *
from commonsdialog import *
Luego, empieza el juego de las líneas de código. Una buena práctica es el reciclaje del código, para los cual definimos las Clases y métodos necesarios.
def cal_suma(mapContext, layer, fieldname):
canv = mapContext.getViewPort()
rect = canv.getEnvelope().getGeometry()
tabla = layer.getFeatureStore()
ents = tabla.getFeatureSet().iterator()
suma = 0
for ent in ents:
g = ent.getDefaultGeometry()
if rect.intersects(g) :
suma += int(ent.get(fieldname))
print "Sumatorio de %s: %s" % (fieldname, suma)
De esta manera obtenemos la vista y su contexto
mapContext = currentView().getMapContext()
Por otro lado así obtenemos la capa activa y solo es cuestión de decir cual campo..
layer = currentLayer()
fieldname = 'Total'
Solo queda llamar a la función para hacer
mapContext = currentView().getMapContext()
…a partir de esto queda poner algunos “adornos” para hacerlo mas interactivo. Para eso usamos el modulo de commonsdialog
El resultado les mostraré en seguida…
En este ejemplo intentaremos calcular un indice de vegetación a partir de una escena Landsat obtenida de http://earthexplorer.usgs.gov/.
Calcularemos el índice de vegetación de diferencia normalizada (NDVI).
Se requiere las bandas que contengan la reflectancia Roja (R ) e Infrarroja cercana (IRC).
La ecuación que usaremos es la siguiente:
\[ NDVI = \frac{IRC-R}{IRC+R} \]
Como resultado obtendremos valores entre -1 y 1
De la API de PyQgis importaremos solo las clases RasterCalculator y RasterCalculatorEntry
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
De la API de PyQgis importaremos solo las clases RasterCalculator y RasterCalculatorEntry
# Función para registrar las capas de entrada
def rasterCalcEntry(EntLst,band,refNam,bNum):
raster = QgsRasterCalculatorEntry()
raster.ref = refNam
raster.raster = band
raster.bandNumber = bNum
EntLst.append(raster)
De la API de PyQgis importaremos solo las clases RasterCalculator y RasterCalculatorEntry
# Definimos el método que calcule el NDVI
def ndvi(R, IRC, output):
Ent = []
rasterCalcEntry(Ent, R, 'R@1', 1)
rasterCalcEntry(Ent, IRC, 'IRC@1', 1)
calc = QgsRasterCalculator(
'(IRC@1-R@1)/(IRC@1+R@1)',
output, 'GTiff', band3.extent(),
band3.width(), band3.height(),Ent)
calc.processCalculation()
Finalmente definimos las entradas y llamamos a las funciones definidas:
# Haremos un listado de capas disponibles
# en el proyecto actual
lyrs = iface.mapCanvas().layers()
# Definirmos la ruta y nombre de salida
output = '../datos/geodata/Landsat/ndvi.tif'
# Llamamos a la función
ndvi(lyrs[1], lyrs[0], output)
Ahora lo veremos en acción…
En este caso no haremos un script de python, sino que emplearemos el lenguaje dentro de una herramienta de un GIS.
La herramienta será la Calculadora de campos de ArcGIS.
El ejemplo fue una práctica dentro del Programa UNIGIS para América Latina.
Consolidaremos los valores de varios campos en uno solo.
Solo escribiremos nuestra función y la ejecutaremos desde la calculadora de campos.
def consolide(*args):
value0 = "-"
for arg in args:
val = str(arg).strip()
if val == "-":
continue
elif val == "":
value0 += " "
else:
value0 += val
return value0.strip("-")
Gabriel V. Gaona
Justus-Liebig Universität Gießen
gabo@gavg712.com
gavg712@gmail.com
gabriel.gaona@umwelt.uni-giessen.de
https://www.researchgate.net/profile/Gabriel_Gaona
Twitter: gavg712