Zonificación Ambiental Energías Renovables, WMS y R parallel
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.
