python

Instalación de librerías NetCDF + Python en un sistema Windows

La instalación adecuada de Python + NetCDF es algo muy extendido y documentado en sistemas Linux, pero no tanto en sistemas Windows. En esta entrada el objetivo es recoger los pasos necesarios para conseguirlo, a modo de referencia rápida con las instrucciones y enlaces oportunos.
Python + netCDF
NetCDF [1] es un formato binario de fichero (.nc) orientado a arrays, que facilita el manejo de varias dimensiones (x, y, z, tiempo…) de forma eficiente y flexible. Por esta razón está muy extendido en la comunidad científica, por ejemplo para guardar datos del medio como temperaturas, corrientes, viento, salinidad, etc. obtenidos mediante sensores o modelado. De forma indisoluble con el formato, existe un conjunto de librerías para NetCDF que facilitan el acceso a sus datos desde varios lenguajes (C, Fortran, Java, etc.) y que son la pasarela con la que trabaja el programador. Python es por su parte el lenguaje de scripting más potente y posee multitud de librerías para ampliar su campo de acción (para el manejo de datos SIG, cálculo intensivo, gráficos, etc.), por lo que juntos forman una buena combinación en el ámbito científico.
Los pasos para configurarlos de forma integrada en una plataforma Windows de 32 bits son:
  • (A) Python 2.7. Si aún no lo tenemos, la última versión disponible para la rama 2.x (la más extendida) es actualmente la 2.7.5, descargable aquí
  • (B) NetCDF4. Las librerías base están escritas en C y para usarlas es necesario descargarlas y compilarlas o, algo más práctico, obtener directamente los binarios. Siguiendo la segunda vía, aquí está el instalador con la última versión para netCDF4. Para una descripción más detallada de estos binarios, ver el enlace [2]
  • (C) Acceso a las librerías NetCDF. Para un enlace posterior a las librerías, es necesario agregar al PATH los siguientes directorios, generados por el instalador: C:\Program Files\netCDF 4.3.0\bin y c:\Program Files\netCDF 4.3.0\deps\shared\bin
  • (D) Numpy. El acceso a los datos de los ficheros NetCDF se realizará de forma efectiva con esta librería para Python, especializada en el manejo de arrays.  La última versión disponible hoy, la 1.7.1, para Python 2.7 está disponible aquí 
  • (E) Interfaz de Python para netcdf4. Ésta es la librería Python que permite leer y escribir en los netcdf desde scripts .py, actuando de pasarela hacia las librerías base previamente configuradas (C). El instalador para la 2.X está en este enlace
  • (F) Validación. Finalmente, para probar que todo está correcto, ejecutaremos un pequeño script de Python, leyendo un netcdf cualquiera (fichero testNC.py). P.ej. estos datos de la NASA muestran las anomalías térmicas sobre la superficie terrestre y es posible descargarlos en formato netCDF desde aquí en un fichero «nmaps.nc». Como nota útil, señalar que para visualizar de forma rápida el netcdf se puede usar un visor como Panoply [3] o un software SIG como gvSIG, que en su versión más reciente incluye soporte a este formato [4]

testNC.py

# -*- coding: utf-8 -*-
import netCDF4 as nc
import numpy as np

'''
     Prueba de acceso a netCDF en Python Win32.
     Abre un fichero .nc y obtiene el valor mínimo de una variable conocida ('TEMPANOMALY')

'''
print 'TEST netCDF en Python'
rutaFichero = "c:/users/admin/nmaps.nc"
fichero = nc.Dataset(rutaFichero)

print "* Variables disponibles en el fichero:"
for v in fichero.variables:
    print v

datos = fichero.variables["TEMPANOMALY"][0]
print "* Mayor anomalia de temperatura negativa : {0} K".format(min(datos))
Nota: como los datos de la variable vienen en un array con «máscara» (para indicar dónde hay ‘huecos’ sin datos), es posible utilizar las clases de numpy específicas y hacer algo como:
minimo = np.ma.MaskedArray.min(fichero.variables["TEMPANOMALY"][:])
maximo = np.ma.MaskedArray.max(fichero.variables["TEMPANOMALY"][:])
Enlaces de interés:

Convertir rasters de un directorio a otro formato con gdal

Script sencillito pero práctico, para convertir todos los ficheros raster de un directorio a otro formato raster (dentro de los soportados por gdal)

Teniendo acceso desde línea de comandos a la utilidad gdal_translate (distribuida p.ej. en FWTools), podemos convertir todos los ficheros .ecw de una carpeta a formato .tif, mediante el siguiente código en Python:

import os, sys
directorioActual = os.getcwd()
directorioSalida = os.path.join(directorioActual, "resultados")

for root, dirs, files in os.walk(directorioActual):
    for fichero in files:
        (nombreFichero, extension) = os.path.splitext(fichero)
        if(extension == ".ecw"):
            comando = "gdal_translate %s %s\\%s" %(fichero, directorioSalida, nombreFichero + ".tif")
            print comando
            os.system(comando)

Copiamos el código a un fichero llamado p.ej. ‘ejecutar.py’ dentro de la carpeta que nos interese, creamos una carpeta dentro para los ficheros transformados (‘resultados’) y ejecutamos… Más fácil imposible.

Se pueden introducir fácilmente otras variantes, p.ej. con otros formatos de raster, viendo las opciones disponibles en la herramienta gdal_translate

EDIT: Editado para añadir mejoras de WordPress al código fuente

Geoproceso interactivo / geoproceso en batch

Geoprocesos ESRIEn ocasiones, con geoprocesos que tienen como parámetro de entrada entidades vectoriales, es interesante que podamos ejecutarlos de dos modos:

(1) de manera interactiva (dibujando nosotros las geometrías) o
(2) suministrando una capa ya existente que tenga las entidades  (p.ej. un shapefile).

Dentro del entorno de geoprocesos de ESRI, esto se corresponde con la definición de un parámetro de tipo FeatureSet (1) y un parámetro como FeatureLayer (2), respectivamente.

[Esto en la versión 9.3 de ArcGIS Desktop viene incorporado ‘de serie’, así que no hay problema. Lo que se comenta a continuación en este post es aplicable solo para la versión 9.2]

Una manera sencilla de que el mismo script pueda trabajar con los dos modos es como se muestra a continuación:

gp = arcgisscripting.create()
gp.overwriteoutput = 1
if gp.ScratchWorkspace is None:
    gp.ScratchWorkspace = gp.GetSystemEnvironment("TEMP")
# Parametros de entrada.........
fs = gp.GetParameter(0)   # (1) FeatureSet [opcional]
fLayer = gp.GetParameter(1)   #(2) FeatureLayer [opcional]
# Detección del MODO DE USO
if (str(fLayer).strip() != ""):
    gp.AddMessage("Trabajando con modo batch...")
else:
    gp.AddMessage("Trabajando con modo interactivo...")
    #Adaptación del FeatureSet a Feature Layer
    tmp = str(int(math.ceil(random.random()*1000000)))
    fLayer = "%s\\%s.shp" %(gp.ScratchWorkspace, tmp)
    fs.Save(fLayer)
# A partir de aquí, el tratamiento puede ser igual
registros = gp.SearchCursor(fLayer)
registro = registros.Next()
while registro:
    # tratamiento del registro [...]
    registro = registros.Next()

En la línea tmp = … también puede usarse gp.CreateScratchName.
Que sea útil.

Descarga de ficheros desde Web Coverage Service (WCS)

Este script sirve para realizar peticiones a un servicio OGC-WCS para, dado un encuadre y una ‘coverage’ conocida, descargar su información a local en ficheros. Está realizado con python y utiliza la librería urllib2. Los ficheros raster generados, en Geotiff o AsciiGrid, pueden luego cargarse en cualquier programa GIS o manipularse con librerías estándar como gdal.

UtilidadesWCS.py

import os
import time
import urllib2

'''
Utilidad de extraccion de servicio WCS
@author: VictorVelarde (victor.velarde@gmail.com)
'''
class ServicioWCS(object):
"""
Utilidad para generar peticiones a un servicio OGC-WCS y guardar los resultados en disco
"""

ANCHURA_PETICION = 10000 #mismas unidades que el sistema de proyeccion
SEGUNDOS_ESPERA_ENTRE_PETICIONES = 5

def __init__(self, url):
''' Constructor con url del servicio, p.ej: http://www.idee.es/wcs/IDEE-WCS-UTM30N/wcsServlet? '''
    self.url = url

def ExtraerFicheros(self, xmin, ymin, xmax, ymax, directorioSalida, formatoWCS):
    if(os.path.isdir(directorioSalida) == False):
        os.mkdir(directorioSalida)

    print 'Inicio de la extraccion'
    i = 1

    for x in range(xmin, xmax, ServicioWCS.ANCHURA_PETICION):
        for y in range(ymin, ymax, ServicioWCS.ANCHURA_PETICION):
            url = ('%s?REQUEST=GetCoverage&SERVICE=WCS&VERSION=1.0.0&FORMAT=%s&COVERAGE=MDT25_peninsula_zip&BBOX=%f,%f,%f,%f&CRS=EPSG:23030&RESX=25&RESY=25' %(self.url, formatoWCS, x, y, x + ServicioWCS.ANCHURA_PETICION,
y + ServicioWCS.ANCHURA_PETICION))
            print url
            peticion = urllib2.Request(url)
            respuesta = urllib2.urlopen(peticion)
            extension = {FormatoWCS.ASCII_GRID : '.asc', FormatoWCS.GEOTIFF: '.tif'}
            fileHandle = open('%s/%03d%s' %(directorioSalida, i, extension[formatoWCS]), 'wb')
            fileHandle.write(respuesta.read())
            fileHandle.close()

            time.sleep(ServicioWCS.SEGUNDOS_ESPERA_ENTRE_PETICIONES)
            i = i + 1

            print 'Fin de la extraccion'

class FormatoWCS:
    ASCII_GRID = 'AsciiGrid'
    GEOTIFF = 'Geotiff'

Y para ejecutarlo:

from UtilidadesWCS import ServicioWCS, FormatoWCS
servicio = ServicioWCS("http://www.idee.es/wcs/IDEE-WCS-UTM30N/wcsServlet")
servicio.ExtraerFicheros(333850.0, 4728100.0, 503850.0, 4888100.0, "extraccion", FormatoWCS.ASCII_GRID)
servicio.ExtraerFicheros(333850.0, 4728100.0, 503850.0, 4888100.0, "extraccion", FormatoWCS.GEOTIFF)

Como alguno me los habéis pedido, aquí subo los ficheros por si alguien quiere probarlos: Megaupload

EDIT: Editado para añadir mejoras de WordPress al código fuente.