Introducción a Python para usos geoespaciales

Roger Veciana i Rovira

Twitter: @rveciana

Blog: http://geoexamples.blogspot.com

Python: Qué es

Lenguaje interpretado
Tipos dinámicos / Tipos fuertes
Orientado a objetos
Multi-propósito
Pensado para ser rápido (de aprendrer, de usar, de entender)
Pensado para crear un código uniforme y limpio

Python: Algunas ventajas

Permite abordar problemas simples de forma simple y problemas complejos
- Programación tipo script
- Programación orientada a objetos con excepciones, etc
Librería estándard muy grande: "batteries included"
Muchas librerías y frameworks que lo extienden
- Creación de entornos web
- Administración de sistemas
- Llibrerías científicas y geogràfiques
- Programación en paralelo
Python une librerías de distintos lenguajes.
Tiene aplicaciones que permiten usarlo de forma interactiva ipython, SAGE

Instrucciones básicas de python

a = 3
b = 4
a + b

Condicionales

if b > a:
    print "b es mayor que a"

Iteraciones

for i in range(10):
    print "Iteración número " + str(i)

NumPy para el tratamiento de matrices

Abrir consola ipython

import numpy
a = numpy.ones((20,10))
a * 23
a + 1
a.sum()
a[3:6,4:8]
a[3:6,4:8] = 2
a[[2,3],[4,4]]

Probar

Otras funciones: max(), mean() (probar a pasar axis=1 de argumento), reshape(), numpy.gradient(a), numpy.arange()
Indexación con otros arrays: a[b], a[b,c]
Operadores lógicos: a>2

GDAL/OGR

Geospatial Data Abstraction Library

Modelo de datos abstracto para leer formatos GIS raster y vectorial
Más de 150 formatos soportados:
- GeoTIFF
- ENVI
- Shapefile
- PostGIS
- GRIB
- NetCDF
- VAISALA IRIS
Usada por muchas aplicaciones y programada en C, tiene bindings en Python
Lectura/Escritura/Cambios de coordenadas/Operaciones GIS

Ejemplo de uso de GDAL

El archivo con los ejemplos es gdal_ejemplos.py

Ejercicios GDAL

Suma de dos arrays
Zonas con más de seis litros en una hora
Zonas donde no llovió a las 00 y si que llovió a las 01
Número de horas de lluvia
Zonas en las que ha llovido todo el dia

Ejemplo de uso de fiona

El archivo con los ejemplos es fiona_ejemplos.py

Ejercicios Fiona

Buscar la distancia entre dos puntos con la función distance()
Si un linestring se crea de la forma: {'type': 'LineString', 'coordinates': ((1, 40),(2,41))}, usar la función crosses e intersection entre dos linestrings
Determinar cuantos rayos han caído en cada comarca (usar comarcas.shp y llamps.shp)

Uso de proyecciones

El archivo con los ejemplos es osr_ejemplos.py

Formato WKT de las proyecciones

'PROJCS["ETRS89 / UTM zone 31N", Nombre de la proyección (arbitratrio)

GEOGCS["ETRS89", Código del Datum
DATUM["European_Terrestrial_Reference_System_1989", Nombre del Datum
SPHEROID["GRS 1980",6378137,298.2572221010002, AUTHORITY["EPSG","7019"]], Definición del elipsoide

TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]], Transformación aproximada del datum a WGS84. 0 indica que el programa asume que son similares
PRIMEM["Greenwich",0],Meridiano de referencia
UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4258"]],Transofmación a radianes

PROJECTION["Transverse_Mercator"],Etiqueta que define la proyección que usará proj4

PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",3],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0], Parámetros de la proyección UTM
UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","25831"]]' Unidades usadas

Geotransform

(258047, 1000, 0.0, 4753640, 0.0, -1000)

Xgeo = gt(0) + Xpixel*gt(1) + Yline*gt(2)
Ygeo = gt(3) + Xpixel*gt(4) + Yline*gt(5)

geotransform = ds.GetGeoTransform()
gdal.ApplyGeoTransform(geotransform, x, y)
gdal.InvGeoTransform(geotransform)

Ejemplo de uso de Mayavi2

El archivo con los ejemplos es mayavi_ejemplos.py