R-Shiny-Data Science

Zonificación Ambiental Energías Renovables, WMS y R parallel

26 Mar 2022 By César Arquero Cabral

Este post está destinado a aquellos que, como en mi caso, os hayais dado de bruces intentando replicar o al menos manejar los nuevos mapas de Zonificación Ambiental para proyectos de energías renovables. En resumen, estos mapas cruzan multitud de fuentes de información para “calificar” el territorio en función de la sensibilidad ambiental. Supuestamente hay unas zonas de valores máximos que deberían ser, bajo los métodos de cálculo, intocables, peeeero…la fuente no es vinculante ni obliga a nada.

¿Cómo se descarga la zonificación?

La respuesta a esta pregunta no la sé. Es posible que se deba a que será un servicio vivo y que puede estar expuesto a modificaciones. Así que, ¿por qué sobrecargar los servidores con miles de descargas que quedarán obsoletas en unos meses? Lo que sí se ofrecen son los servicios WMS para cargar en vuestros proyectos:

Si no es suficiente el WMS para ti, entonces tendrás que trabajar un poco. Vamos allá.

Capturar el servicio WMS (método QGIS a lo bruto)

QGIS puede ser utilizado de forma un poco salvaje para descargar un servicio WMS completo. Claro, son imágenes, pero si calculáis bien la resolución y dimensiones totales y de las partes del raster podréis conseguir un buen resultado. Quizá el proceso de cálculo de teselas, generación del atlas y exportación podrían ser explicadas en otro artículo. De momento lo indico en un pequeño resumen:

Diseñar un mapa sencillo en el atlas. Si es demasiado grande (por ejemplo la extensión de España) podréis ser bloqueados por el servidor o bloqueados por vuestos modulos de RAM. En mi caso, el diseño de hoja es de 1x1m. Not bad.
Calcular la resolución y escala a la que la imagen empieza a “mostrar” el tamaño de pixel. Lo podéis ver también en la información del servicio.
Crear un grid con vuestra herramienta favorita a partir de la escala total y escala a la que irán las hojas tendréis más o menos teselas. Estas serán las teselas a imprimir.
El grid servirá de atlas. Así que, una vez tengáis todo listo, activáis el atlas y vais exportando a .tif o a otro formato pero con información georreferenciada. Haced esto por la noche que los servidores estarán más ligeros y no molestaréis a nadie.
Construir un raster virtual a partir de todas esas hojas generadas y listo. Ya tenéis el servicio descargado.

Aquí una pequeña vista para que os hagáis una idea.

Un momento… ¡¡Malditos raster multibanda, solo quiero los valores!!

Este es el meollo del artículo. Malditos colores y bandas, si tan siquiera pudiésemos tener simples valores numéricos. ¿Cómo se hace eso? Volviendo a trabajar.

La transformación de los rasters puede ser realizada en R, como todo, como siempre. En este caso os expongo el código completo junto a un pequeño resumen de lo que se ejecuta en él. El código es suficientemente explicativo, aunque cabe destacar que hay partes que se ejecutan con procesos en paralelo y están configuradas para ejecutarse en entorno linux. Por lo demás, nada de qué preocuparse.

Procesado de rasters en paralelo

El resumen del proceso:

Establecer los parámetros: directorios donde se alojan los archivos raster georreferenciados, las carpetas de salida para los archivos transformados en valores y nombre del raster virtual a generar con ellos. Calcular las combinaciones posibles de colores. Esto es necesario para establecer el gradiente correcto de valores. En este proceso básicamente se leerán los rasters en paralelo y se extractarán todos los colores diferentes en cada uno. Los valores se convertirán en formato hex.
Ordenar colores. Es uno de los principales hándicaps de haber descargado el wms a lo bruto. A partir del visor web se ha comprobado que los valores van de 0 a 10, siendo 10 el máximo valor ambiental. En base a esto, se asignarán los valores a los colores. En este caso, el naranja más fuerte será el 10.
Transformar los colores de los rasters en valores 0 a 10. Esto se realiza también en paralelo para agilizar el proceso. Esto permite (1) aligerar los archivos iniciales al pasar de un raster multibanda a un raster singleband y (2) poder fliparse.
Crear un raster virtual con los nuevos archivos. Una función muy sencilla de utilizar está en la librería landtools, todavía en proceso de pulido.

library(raster)
library(ggplot2)
library(parallel)

# Parámetros
# Realizar cambios en función de la fuente y localización de los archivos raster
# generados a partir del proceso explicado para QGIS
source <- '../01_Fuentes/02_GIS/Zonificación/EAEF/E25000_CAN/'
r <- stack('../01_Fuentes/02_GIS/Zonificación/EAEF/EAEF_E25000_CAN.vrt')
outdir <- '../01_Fuentes/02_GIS/Zonificación/EAEF/E25000_CAN_val/'
outname <- 'EAEF.tif'

# Salida para el vrt
out <- '../01_Fuentes/02_GIS/Zonificación/EAEF/'
vrtname <- 'EAEF_E25000_CAN_val.vrt'

#-------------------------------------------------------------------------------
# 1. Calcular todas las combinaciones posibles de colores RGB
#-------------------------------------------------------------------------------
files <- list.files(source, '.tif$', full.names = TRUE)

#-----
# Proceso en paralelo con barra de progreso
library(doSNOW)
library(tcltk)

cl <- makeSOCKcluster(detectCores())
registerDoSNOW(cl)

pb <- txtProgressBar(max=100, style=3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress=progress)
vals <- foreach(i=seq_along(files), .options.snow=opts) %dopar% {
        unique(paste(raster::values(raster::raster(files[[i]], 1)),
                     raster::values(raster::raster(files[[i]], 2)),
                     raster::values(raster::raster(files[[i]], 3))))
}
close(pb)
stopCluster(cl)
#-----
RGBcomb <- unique(unlist(vals))

# Agrupar los valores
numbs <- str_split(RGBcomb, ' ')

# Crear lista de valores en formato hex
hexs <- list()
for(i in seq_along(numbs)){
        x <- as.numeric(numbs[[i]])
        hexs[[i]] <- rgb(x[1], x[2], x[3], maxColorValue = 255)
}

#-------------------------------------------------------------------------------
# 2. ORDENAR COLORES, 
# https://stackoverflow.com/questions/61193516/how-to-sort-colours-in-r
# Por suerte solo hace falta ordenarlos por intensidad
#-------------------------------------------------------------------------------
library(TSP)
rgbs <- col2rgb(hexs)
lab <- convertColor(t(rgbs), 'sRGB', 'Lab')
oc2 <- hexs[order(lab[, 'L'])]
qplot(x = seq_along(hexs), y = 1, fill = I(oc2), geom = 'col', width = 1) +
    theme_void()

# color alto = 10, color bajo = 0 (sin valor)
# el resto de valores se dividen por igual  
maxhex <- oc2[[1]]
minhex <- oc2[[length(oc2)]]
maxval <- 10
minval <- 0
interhex <- length(oc2)-2
interdiv <- ceiling(interhex/9)
vals <- rep(9:1, times = interhex, length.out = interhex, each = interdiv)
vals <- c(maxval, vals, minval)

# Compound list of values (index)
hex <- unlist(oc2)
val <- vals

#-------------------------------------------------------------------------------
# 3. Transformar todos los rasters de colores en raster de valores
#-------------------------------------------------------------------------------
cl <- parallel::makeCluster(detectCores(), type="FORK")
doParallel::registerDoParallel(cl)
foreach(i=seq_along(files)) %dopar% {
        r <- stack(files[i])
        rasterhex <- rgb(values(r), maxColorValue = 255)
        
        # map values
        nval <- as.numeric(plyr::mapvalues(rasterhex, from=hex, to=val))
        
        # Create a raster replacing hex values for vals calculated
        newr <- raster()
        res(newr) <- res(r)
        dim(newr) <- dim(r)
        crs(newr) <- crs(r)
        extent(newr) <- extent(r)
        
        # Change raster values
        newr[] <- nval
        
        # Export raster
        # fname <- str_replace(basename(files[i]), '.tif$', '')
        fname <- basename(files[i])
        writeRaster(newr, paste(outdir, fname), datatype="FLT4S", overwrite=T)
}
stopCluster(cl)

#-------------------------------------------------------------------------------
# 4. Construir un nuevo raster virtual con los nuevos archivos raster de valores
#-------------------------------------------------------------------------------
library(landtools)
library(rgdal)
library(gdalUtils)
landtools::dir2vrt(outdir, out, vrtname = vrtname, pattern = '.tif$')

Aquí podéis ver cómo la ejecución de los procesos en paralelo os permite generar archivos raster modificados a toda leche. Os lo recomiendo mucho para ejecutar cualquier tipo de tareas sobre ellos (crops, masks, reclass, etc). Con esto estarían listos los rasters donde como mínimo podréis ver las zonas con valores máximos y qué proyectos, presentes o futuros, están actualmente sobre ellas.