Archivo de la categoría: Tutoriales R

¿Cómo descargar información de GBIF usando R?

Muchas investigaciones se realizan en identificar cuál es el impacto del cambio climático sobre ciertas especies tanto animales como vegetales, dentro de los métodos para dichas investigaciones se encuentran MaxEnt y Random Forest, para hacer tales modelaciones se hace necesario contar con presencias de la especie a modelar. GBIF (Global Biodiversity Information Facility) es una plataforma que cuenta con cientos de miles de registros de presencias tanto de flora como de fauna, estos registros en su mayoría cuentan con las coordenadas geográficas que pueden ser de uso para realizar el tipo de investigaciones mencionado anteriormente.

Para hacer la descarga de tales registros se puede hacer de varias maneras, una de ellas es haciendo el respectivo registro en la página web (requiere usuario y contraseña), y luego introduciendo el nombre científico de la especie nos arroja el resultado de todos los registros, sin embargo, en muchas ocasiones el hacer la descarga se demora desde pocos minutos hasta varias horas (dependiendo del número de registros). Otro camino para descargar la información es hacer uso de la libreria RGBIF del lenguaje de programación R, aquí es mucho más rápido hacer la descarga de información, sumado a que ya tendríamos la tabla cargada en R para hacer depuración de datos, análisis exploratorio espacial, entre otros.

En este tutorial entonces se explicara como hacer la descarga de registros de la especie Musa paradisiaca haciendo uso de R, y luego como obtener unicamente los registros localizados en el país de Colombia.

Se sugiere tener un nivel básico en el manejo de R para realizar el ejercicio, sin embargo, se espera ser lo más claro posible para que pueda ser replicado sin mayor problema. Se require tener instalado el software R Studio.

Cómo primer paso se descarga la libreria de RGBIF, tambien se instalara la libreria Tidyverse (para hacer manejo de la tabla de una manera más facil).

install.packages('rgbif')
install.packages('tidyverse')
library(rgbif)
library(tidyverse)

Ahora se hará la respectiva descarga de los datos para la especie Musa Paradisiaca (conocida popularmente como banano).

occ <- occ_data(scientificName = 'Musa paradisiaca L.', 
limit = 200000,
hasCoordinate = TRUE,
hasGeospatialIssue = FALSE)

Ahora conoceremos los nombres de las columnas. Las tablas descargadas de RGBIF siguen la estructura de darwin core (conocer más de este formato aquí).

colnames(occ$data)

Una de las columnas incluye el nombre del país (country), entonces, con base en esta columna haremos el respectivo filtro para obtener los registros localizados unicamente en Colombia.

occ_col <- filter(occ$data, country == 'Colombia')

Con esto entonces, el objeto occ_col contiene unicamente las presencias de la especie Musa Paradisiaca para el país de Colombia. Ahora hacemos la descarga del shapefile de los departamentos de Colombia, y se procede a la realización de un mapa sencillo, para visualizar en donde se encuentran los puntos de Musa paradisiaca en Colombia según esta base de datos de GBIF.

shp <- raster::getData('GADM', country = 'COL', level = 1)

gg <- ggplot() +
geom_point(data = occ_col, aes(x = decimalLongitude, y = decimalLatitude), color = 'black') +
geom_polygon(data = shp, aes(x = long, y = lat, group = group), color = 'grey', fill = NA) +
coord_fixed(ylim = c(-5, 12.5), xlim = c(-80, -67)) +
xlab('Longitud') +
ylab('Latitud')

ggsave(plot = gg, filename = 'myMap.png', units = 'cm', width = 12, height = 16, dpi = 300)

Ahora finalizamos haciendo la escritura de la tabla en nuestro espacio de trabajo:

write.csv(occ_col, 'occ_musa_paradisiaca.csv', row.names = F)

El mapa que generamos con el objeto gg debería ser similar al siguiente:

Creando un mapa en R

Hola a todos!

En este tutorial realizaremos un mapa de una variable continua en R, en el cuál utilizaremos archivos espaciales tipo raster, y shapefile. Este mapa representará la ubicación geográfica de estaciones meteorológicas de precipitación, así como tambien la temperatura promedio para el departamento del Valle del Cauca. En caso que quiera replicar el ejercicio puede descargar los datos aquí.

Aquí se hacen uso de algunos comandos como filter, select, extract, which, getData.

El primer paso a realizar es cargar las librerias, aquí usamos el comando p_load de la libreria pacman para cargar varias librerias en una sola línea, además, en caso que la libreria no encuentre instalada, ésta se instalará automaticamente.

require(pacman)
pacman::p_load(raster, rgdal, tidyverse, rgeos, gtools, stringr, foreign)

El segundo paso es cargar la tabla con las coordenadas geográficas de las estaciones de Worldclim.

wcs <- read.dbf('_dbf/wc_tean_stations.dbf')

Como tercer paso esta descargar el límite administrativo del Colombia a nivel administrativo 1 y 2.

col <- raster::getData('GADM', country = 'COL', level = 1)
col2 <- raster::getData('GADM', country = 'COL', level = 2)
tmn <- raster::getData('worldclim', var = 'tmean', res = 0.5, lon = coordinates(vll)[[1]], lat = coordinates(vll)[[2]])

El cuarto paso es realizar la extracción por máscara unicamente para el Valle del Cauca; pues los datos descargados corresponden a una ventana que contiene todo Colombia y algunos países centroamericanos.

vll <- col[col@data$NAME_1 %in% 'Valle del Cauca',]
vl2 <- co2[co2@data$NAME_1 %in% 'Valle del Cauca',]
tmn <- raster::crop(tmn, vll) %>% raster::mask(vll)

Ahora si realizamos un plot sencillo de los datos raster podemos observar que hay una grande zona en la que no hay datos, aquí se podría hacer una exclusión de estos datos NA del raster, siguiendo este proceso:

plot(tmn)
xNA <- is.na(tmn)
colNotNA <- which(colSums(xNA) != nrow(x))
rowNotNA <- which(rowSums(xNA) != ncol(x))
croppedExtent <- extent(tmn,
r1 = rowNotNA[1],
r2 = rowNotNA[length(rowNotNA)],
c1 = colNotNA[1],
c2 = colNotNA[length(colNotNA)])
tmncut <- raster::crop(tmn, croppedExtent)

Ahora bien, procedemos a realizar el calculo del promedio de la temperatura de los 12 meses del año y a su vez se divide el raster entre 10 (pues los datos en crudo vienen multiplicados por 10), esto mediante las siguientes  dos líneas de comando.

tmn.avg <- mean(tmncut)
tmn.avg <- tmncut / 10

El siguiente paso es hacer el filtro de las estaciones climáticas únicamente para Colombia, a su vez solo escogemos las columnas country, name, long y lat.

wcs <- wcs %>% 
as.tibble() %>%
dplyr::filter(COUNTRY == 'COLOMBIA') %>%
dplyr::select(COUNTRY, NAME, LONG, LAT)

Ahora el siguiente paso es realizar la extracción de los nombres de los departamentos para las estaciones, para eso hacemos uso del shape de nivel administrativo 2 así como tambien de la tabla de las estaciones meteorológicas

wcs <- raster::extract(co2, wcs[,c('LONG', 'LAT')]) %>% dplyr::select(NAME_1, NAME_2) %>% cbind(wcs, .) 
wcs <- wcs %>% dplyr::filter(NAME_1 == 'Valle del Cauca')

Ahora sí, procedemos a realizar el mapa haciendo uso de funciones de la libreria ggplot.

gg <- rasterVis::gplot(tmn.avg) +
geom_tile(aes(fill = value)) +
scale_fill_gradientn(colours = RColorBrewer::brewer.pal(n = 8, name = "YlOrRd"), na.value = 'white') +
labs(fill = '') +
geom_polygon(data = vl2, aes(x=long, y = lat, group = group), color = 'grey', fill='NA') +
geom_polygon(data = col, aes(x=long, y = lat, group = group), color = 'grey', fill='NA') +
coord_equal(xlim = c(-77.6, extent(vl2)[2]), ylim = extent(vl2)[3:4]) +
xlab('Longitud') +
ylab('Latitud') +
theme_bw() +
theme(panel.background = element_rect(fill = 'white', colour = 'black'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = 'top')

gg <- gg +
geom_point(data = wcs, aes(x = LONG, y = LAT), col = 'blue')

Sí queremos visualizar el mapa escribimos

gg

Por último, procedemos a hacer el guardado de nuestro mapa:

ggsave(plot = gg, filename = 'myMap.png', unit = 'in', width = 9, height = 7, dpi = 300)

Descarga y tratamiento de datos climáticos a nivel mensual (CRU) usando R

En este tutorial trabajaremos con datos climáticos descargados de la pagina web de CRU (Climatic Researc Unit), se realizará la descarga, y tratamiento de los datos espaciales con el fin de realizar gráficas de tendencia a lo largo de 35 años (1980-2015); todo lo anterior lo realizaremos en el lenguaje de programación R. Las variables disponibles para descarga son las siguientes:

– cld: cobertura de nubes (%)
– dtr: rango temperatura diurno (grados centigrados)
– fdf: frost day frequency (días)
– pet: evapotranspiración potencial (mm/day)
– pre: precipitation (mm por mes)
– rhm: humedad relativa (porcentaje)
– ssh: sunshine duración (hours)
– tmp: temperatura promedio (grados centigrados)
– tmn: promedio de temperatura minima (grados centigrados)
– tmx: promedio de temperatura máxima (grados centigrados)
– presión de vapor (hectopascal hPa)
– wet: frecuencia de días húmedos (días)
– wnd: velocidad del viento (metros por segundo)

En este ejercicio usaremos las variables de precipitación, temperatura máxima, promedio y mínima. Para la realización de este ejercicio utilizaremos información que descargaremos durante el mismo tutorial. Se recomienda tener un nivel básico de programación para poder entender este tutorial, sin embargo, se intentará explicar las lineas a utilizar.

A continuación, listaremos los pasos a realizar explicando cada uno de ellos.

1. Cómo primer paso esta la instalación de los paquetes a utilizar en el sitio web.

# Load libraries
if(!require(raster)){install.packages('raster'); library(raster)} else {library(raster)}
if(!require(raster)){install.packages('rgdal'); library(raster)} else {library(rgdal)}
if(!require(tidyverse)){install.packages('tidyverse'); library(tidyverse)} else {library(tidyverse)}
if(!require(utils)){install.packages('utils'); library(utils)} else {library(utils)}
if(!require(gtools)){install.packages('gtools'); library(gtools)} else {library(gtools)}
if(!require(velox)){install.packages('velox'); library(velox)} else {library(velox)}
if(!require(stringr)){install.packages('stringr'); library(stringr)} else {library(stringr)}

2. Luego procedemos a configurar las variables y los años a descargar, para ello creamos los objetos vars, yr.start y yr.end; además establecemos la ruta base.

path_base <- 'https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.01/cruts.1709081022.v4.01/'
vars <- c('pre', 'tmn', 'tmx', 'tmp')
yr.start <- c(1981, 1991, 2001, 2011)
yr.end <- c(1990, 2000, 2010, 2016)

3. Ahora procedemos a crear dos ciclos, haciendo uso del comando lapply, para descargar tanto las 4 variables como los distintos periodos de tiempo, de manera automática.

lapply(1:length(vars), function(v){
print(vars[v])
lapply(1:length(yr.start), function(x){
print(yr.start[x])
download.file(url = paste0(path_base, vars[v], '/cru_ts4.01.', yr.start[x], '.', yr.end[x], '.', vars[v], '.dat.nc.gz'),
destfile = paste0('_raster/_gz/', vars[v], '_', yr.start[x], '_', yr.end[x], '.nc.gz'),
mode = 'wb')
})
})

4. Como siguiente paso procedemos a descomprimir los archivo .gz obtenidos, para ello hacemos uso de las siguientes lineas. En el parametro destname establecemos la ruta donde queremos guardar nuestro conjunto de datos.

Todos los archivos resultantes los guardaremos en la carpeta nc (parámetro destname de la línea en la que usamos el comando gunzip).

gzfiles <- list.files('../_download/_raster/_gz', full.names = T, pattern = '.gz$')
nms <- basename(gzfiles) %>% gsub('.gz', '', .)
lapply(1:length(gzfiles), function(x){
print(gzfiles[x])
R.utils::gunzip(filename = gzfiles[x], destname = paste0('../_download/_raster/_nc/', nms[x]))
})

5. Ahora procedemos a listar los objetos que hemos descomprimido y luego las decadas de dichos archivos, estos elementos serán los parámetros a utilizar en la función que hará la extracción por máscara para Colombia de los archiovs, y luego la escritura de los mismos en una carpeta de nuestro espacio de trabajo.

ncfiles <- list.files('../_download/_raster/_nc', full.names = T, pattern = '.nc$')
decs <- basename(ncfiles) %>% str_sub(., start = 5, end = 13) %>% mixedsort() %>% unique()
col <- getData('GADM', country = 'COL', level = 0)
extNC <- function(vr, dec){
ncfile <- grep(vr, ncfiles, value = T) %>% grep(dec, ., value = T)
stk <- brick(ncfile)
stk.col <- raster::crop(stk, col) %>% raster::mask(., col)
nms <- names(stk.col) %>% gsub('\\.', '_', .) %>% str_sub(string = ., start = 2, end = 8)
lyrs <- unstack(stk.col)
Map('writeRaster', x = lyrs, filename = paste0('../_download/_raster/_tif/', vr, '_', nms, '.tif'))
print('Done!')
}

Aquí la aplicación de la función anteriormente creada. Este ciclo intera primero en la variable y luego en la decena, así primero extraera todos los períodos de tiempo y luego pasará a la siguiente variable, hasta completar las cuatro variables descargadas.

lapply(1:length(vars), function(v){
print(vars[v])
lapply(1:length(decs), function(d){
print(decs[d])
extNC(vr = vars[v], dec = decs[d])
})
})

6. En la función anterior logramos obtener, entonces, los archivos de precipitación, temperatura máxima, media y mínima para Colombia. Ahora plotearemos un archivo de ellos, y después convertiremos los archivos raster en tablas para de esta manera hacer gráficas de tendencia para las variables mencionadas anteriormente.

Lo primero que haremos en este paso será listar los archivos resultantes del pasos anterior, y luego a partir de un for crearemos 4 stacks, uno para cada variable, tal como se muestra a continuación:

fls <- list.files('../_download/_raster/_tif', full.names = T, pattern = '.tif$') %>% mixedsort()
for(i in 1:length(vars)){
eval(parse(text = paste0(vars[i], '<- grep(vars[', i, '], fls, value = T) %>% stack()')))
}
plot(pre[[1]])

7. Lo siguiente a realizar corresponde a la conversión de los archivos raster a tabla, para poder así crear las gráficas de tendencias. Aquí crearemos la función que realiza dicha conversión.

Lo primero que se realiza en la función es la conversión del stack a un objeto velox, esto pues es mucho más rapido las funciones de dicha libreria, luego extraemos las coordenadas de los raster, y luego con SpatialPoints convertimos las coordeanas en un shapefile, el siguiente paso es la extracción de los valores para todos y cada uno de los raster, luego y quizá estas son las lineas más complejas de entender, es la conversión de la matrix resultante en un data frame, aquí se realiza la asignación de un ID para cada píxel, así como también la creación de columnas que representan el año y el mes.

rstToTbl <- function(rst){
print(rst)
#rst <- pre
vlx <- velox(rst)
crd <- vlx$getCoordinates()
spn <- SpatialPoints(coords = crd)
vls <- vlx$extract_points(sp = spn)
rsl <- vls %>%
tbl_df() %>%
.[complete.cases(.),] %>%
mutate(ID = 1:nrow(.)) %>%
setNames(c(names(rst), 'ID')) %>%
gather(var, value, -ID) %>%
mutate(year = str_sub(var, 5, 8),
mnth = str_sub(var, 10, 11))
print('Done')
return(rsl)
}

Ahora aplicamos la función a las variables, para ello primero creamos una lista a partir de los objetos: pre, tmx, tmp, y tmn. Y seguido de ello con el ciclco del lapply aplicamos la función a la lista, el resultado será el objeto fnl (de final) que contendrá los cuatro tablas con la información climática de Colombia.

stks <- list(pre, tmx, tmp, tmn)
fnl <- lapply(1:length(stks), function(k) rstToTbl(rst = stks[[k]]))

8. Ahora haremos uso de funciones de la libreria ggplot para crear una gráfica para un píxel en particular. Sí Ud. lo desea puede crear un ciclo (lapply o for, por ejemplo) para crear una gráfica para cada píxel.

En este paso lo primero que haremos será un resumen anual de la precipitación, es decir, haremos la suma de la precipitación para los 12 meses de cada año, esto para un píxel en particular (ID = 300); luego crearemos el gráfico y estimaremos una gráfica de tendencia, vale aclarar que existen muchas maneras de realizarlo, este es solamente uno de ellas.

x <- fnl[[1]] %>%
filter(ID == 300)
x <- x %>%
group_by(year) %>%
summarize(sum = sum(value)) %>%
mutate(year = as.numeric(year))
gg <- ggplot(data = x, aes(x = year, y = sum)) +
geom_line() +
geom_smooth(method = 'lm') +
xlab('year') +
ylab('Prec (mm)') +
theme_bw()
ggsave(plot = gg, filename = 'prec_line.png', units = 'cm', width = 15, height = 13, dpi = 150)

El resultado debería ser algo como esto:

Sí has llegado hasta gracias por leer el tutorial, y por favor, no olvides compartirlo para que así este sea difundido de la mejor manera.

Conversión de Archivos KML a Shapefile mediante procesamiento en Python y R

En este tutorial explicare el proceso de importar archivos KML a shapefile mediante dos códigos, el primero escrito en el lenguaje de programación Python, y el segundo en el lenguaje R; en el primer proceso se exportan cada archivo Kml a una geodatabase, y luego en R importamos los archivos de la geodatadase y los exportamos como un archivo shapefile. Para llevar a cabo el presente tutorial es necesario tener instalado cualquier versión de ArcMAP superior a 10, y RStudio versión 3 o superior.

La ventaja que nos brinda hacer el proceso mediante línea de código es que si tenemos, por ejemplo, 100 archivos KML no sería necesario hacer el proceso 100 veces con la función ‘Import KML to Layer’ de ArcGIS, sino de una manera mucho más eficiente con un ciclo en ambos lenguajes de programación.

Los archivos para realizar el presente tutorial están disponibles en este enlace web. Estos archivos corresponden a polígonos dibujados aleatoriamente en Google Earth, no corresponden a ningún tipo de cobertura, uso de la tierra, o elemento geográfico; fueron realizados en Google Earth solamente para fines de hacer este tutorial.

Aquí los pasos a seguir:

Sección Python

1. En un compilador de texto, escribimos las siguientes líneas de código, el compilador que yo utilizare será notepad (clic aquí para descargar). Las primeras líneas son las siguientes:

import arcpy, os
arcpy.env.workspace = r'path\_kmz'
outLocation = r'path\_gdb'
MasterGDB = 'AllKMLLayers.gdb'
MasterGDBLocation = os.path.join(outLocation, MasterGDB)

En la segunda línea indicamos cuál será nuestro espacio de trabajo, en este espacio de trabajo es donde se encuentra el listado de archivos KMZ. En la tercer linea va la ruta donde quedaran guardadas las geodatabases. Y por último, en la cuarta línea escribimos el nombre de cómo queremos que llame la geodatabase donde guardaremos los layers.
2.      Ahora escribiremos estas líneas, aquí lo que haremos es hacer un for para leer los archivos .KM que contiene la carpeta de nuestro espacio de trabajo, y establecemos el ambiente de trabajo donde se encontraran nuestros archivos de salida.
print 'Master'

for kmz in arcpy.ListFiles('*.KM*'):
print ("CONVERTING: {0}".format(os.path.join(arcpy.env.workspace, kmz)))
arcpy.KMLToLayer_conversion(kmz, outLocation)

arcpy.env.workspace = outLocation
wks = arcpy.ListWorkspaces('*', 'FileGDB')
wks.remove(MasterGDBLocation)
1.   3.   El siguiente paso será escribir el último grupo de líneas, en donde se ejecuta la conversión de los archivos kmz a geodatabase mediante un for, este permite iterar en caso dado que se tengan muchos archivos a convertir

for fdgb in wks:

arcpy.env.workspace = fdgb
featureClasses = arcpy.ListFeatureClasses('*', '', 'Placemarks')

for fc in featureClasses:
print ("COPYING: {0} FROM: {1}".format(fc, fgdb))
fcCopy = os.path.join(fgdb, 'Placemarks', fc)
arcpy.FeatureClassToFeatureClass_conversion(fcCopy, MasterGDBLocation, fgdb[fgdb.rfind(os.sep)+1:-4])


# Clean up
del kmz, wks, fc, featureClasses, fgdb

Estas líneas de código fueron escritas a partir de la ayuda que ofrece ArcMap (sí desea visitar la página web aquí el enlace).
Sección R
Ahora pasaremos a abrir RStudio y escribiremos las siguientes líneas.
1.   1. El primer pasos será escribir, de igual manera como en Python, la carpeta de trabajo en donde se encuentran las geodatabases creadas mediante el punto anterior. Antes de ello es importante cargar las librerías a utilizar, en caso dado que usted no tenga las librerías las puede instalar mediante el comando install.packages (‘nameLibrary’).
require(raster)
require(rgdal)
require(plotKML)
require(maptools)
require(XML)

files <- list.files('_data/_gdb', full.names = T, pattern = '.gdb')
list <- lapply(files, FUN = ogrListLayers)

En la línea del lapply aplicamos la función ogrListLayers a las geodatabases, esto con la finalidad de conocer cuál es el nombre de cada shapefile dentro de la geodatabase.
2. Luego creamos un objetos tipo lista que es donde guardaremos los archivos con los shapefiles
poly  <- list()
1.  
      3. Seguido de ello pasamos a escribir un ciclo for, mediante el cuál leemos los archivos shapefiles dentro de la geodatabase, y escribimos estos archivos en una carpeta distinta a la que contiene las geodatabases. Dentro del for hay un condicional, este sirve para conocer si el archivo que vamos a escribir ya existe dentro de nuestro espacio de trabajo, en caso dado que ya exista el archivo no lo escribe de nuevo. 

for(i in 1:length(files)){

poly[[i]] <- readOGR(dsn = files[[i]], layer = as.character(list[[i]][1]))

if(!file.exists(paste0('/_data/_shp/', gsub(' ', '', as.character(poly[[i]]@data[1,1])), '.shp'))){

writeOGR(obj = poly[[i]], dsn = '_data/_shp', layer = gsub(' ', '', as.character(poly[[i]]@data[1,1])), driver = 'ESRI Shapefile')
print('Written Shapefiles')

}

}

Los códigos los puedes encontrar, además, en mi cuenta de Github, disponible en este enlace.

Realizar extracción por máscara en R haciendo uso de paralelización para optimizar el tiempo en el procesamiento

En este tutorial realizaremos una extracción por máscara de 48 archivos raster, que representan la precipitación, temperatura máxima, media y mínima, para 3 distintos departamentos de Colombia; la idea final es obtener los 48 raster para cada departamento, es decir que se obtendrían al final 144 raster. Para ello haremos uso de paralelización, lo que permite ejecutar el mismo proceso para los tres distintos departamentos en núcleos distintos, es decir, a la misma vez que se cortan los datos para Antioquia se están cortando los datos para Cundinamarca y Valle del Cauca, pues se hará uso de tres núcleos (si su computador tiene menos de cuatro núcleos no podrá replicar este ejercicio tal cual).

Los tres departamentos a los cuales les haremos la extracción por máscara son Antioquia, Valle del Cauca y Cundinamarca, para ello tendremos el shapefile de división político administrativa de Colombia (Fuente: SIGOT); además de ello, los raster climáticos para toda Colombia. El ejercicio incluye la descarga de los datos para que replique el ejercicio (URL).

1. Cargamos las librerías necesarias, en caso dado que usted no tenga la librería esta se descargara automáticamente.

if(!require(raster)){install.packages('raster'); library(raster)} else {library(raster)}
if(!require(rgdal)){install.packages('rgdal'); library(rgdal)} else {library(rgdal)}
if(!require(spdplyr)){install.packages('spdplyr'); library(spdplyr)} else {library(spdplyr)}
if(!require(tidyverse)){install.packages('tidyverse'); library(tidyverse)} else {library(tidyverse)}
if(!require(foreach)){install.packages('foreach'); library(foreach)} else {library(foreach)}
if(!require(parallel)){install.packages('parallel'); library(parallel)} else {library(parallel)}
if(!require(doSNOW)){install.packages('doSNOW'); library(doSNOW)} else {library(doSNOW)}
if(!require(gridExtra)){install.packages('gridExtra'); library(gridExtra)} else {library(gridExtra)}
if(!require(rasterVis)){install.packages('rasterVis'); library(rasterVis)} else {library(rasterVis)}

2. Descargamos y declaramos los archivos a utilizar en este tutorial, tanto los archivos raster como el shapefile de departamentos de Colombia

url    <- 'https://zdl9fa.bn1303.livefilestore.com/y4mV948RDpIKmA015QmtTAON46y0O8vW0kF7B3JQ6xeP006WJDddpGmwQ1wpRS9G3f36LGwVNy-EbwxqA9HPNUX6msaosaUNe3RJf8WynIjTbMohC490Mer7l5lrpEMTdVPR1-XWmPHv9CUaPxLFnqqu1qLIo5a36EjDeR4Ry0QmCsvt4G9TXyoLVoMkWt54VOUpAH4VIsLSuS5yk44wVI8Iw/_data.zip'
mydir  <- 'D:/parallel'
if(!file.exists(mydir)){dir.create(mydir)}
temp <- tempfile(tmpdir = mydir, fileext = '.zip')
download.file(url, temp)
unzip(temp, exdir = mydir)
unlink(temp) # Eliminar el archivo .zip

adm <- shapefile(paste0(mydir, '/_data/_shp/LimiteDptal_geo_ok.shp'))
lyrs <- list.files(paste0(mydir, '/_data/_raster/_co'), full.names = T, pattern = '.asc$') %>%
lapply(FUN = raster) %>%
stack()

3.  Del shape de Colombia extraemos tres objetos distintos que sean cada departamento, y luego los guardamos en un objeto tipo lista. 

ant   <- filter(adm, NOMBRE_DPT == 'ANTIOQUIA')
val <- filter(adm, NOMBRE_DPT == 'VALLE DEL CAUCA')
cun <- filter(adm, NOMBRE_DPT == 'CUNDINAMARCA')
shps <- list(ant, val, cun)

4. Creamos las carpetas donde se guardaran los archivos cortados para cada departamento, en mi caso creare ‘_ant’ para Antioquía, ‘_val’ para Valle del Cauca, ‘_cun’ para Cundinamarca.

dirsOut <- paste0(path, '/', c('_ant', '_val', '_cun'))
if(!file.exists(dirsOut)){lapply(dirsOut, FUN = dir.create)}

5.  Indicamos la cantidad de núcleos a utilizar. Se recomienda utilizar un núcleo menos que el total de núcleos que tenga su computador.

nCores   <- detectCores(all.tests = FALSE, logical = TRUE)
cl <- makeCluster(nCores - 1) #Número de nucleos a utilizar
registerDoSNOW(cl)

6. Configuramos la barra de progreso para así conocer que tanto hace falta para que se termine el proceso de extracción por máscara.

length_run <- length(shps)
pb <- txtProgressBar(max = length_run, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress=progress)


7. Creación de la función de extracción por máscara. Primero usamos la función, crop y mask para cortar por departamento y luego realizamos un unstack para tener los archivos raster como lista, y poderlos así escribir mediante un for sencillo.

extractMask <- function(lyrs, shps, dirsOut){

lyrs_cut <- crop(lyrs, shps) %>%
mask(shps) %>%
unstack()

for(j in 1:length(lyrs_cut)){

writeRaster(lyrs_cut[[j]], paste0(dirsOut, '/', names(lyrs_cut[[j]]), '.asc'))

}

return(lyrs_cut)

}

8. Ejecutamos el foreach para correr el proceso paralelamente en los tres núcleos destinados para ello, se recomienda no tener abierto muchos programas en el computador para que así no se ponga lento el PC.

lyrs_cut_ok <- foreach(i = 1:length(shps), .packages = c('raster', 'rgdal', 'dplyr'), .options.snow=opts) %dopar% {

extractMask(lyrs, shps[[i]], dirsOut[i])

}

9. Creamos una función sencilla para plotear los archivos resultantes

plotGrid <- function(lyr){

gg <- gplot(lyr) +
geom_tile(aes(fill = value)) +
coord_equal() +
labs(fill = "Precipitación \n Enero (mm)") +
xlab('Lon') +
ylab('Lat')

}

plotAnt <- plotGrid(lyr = lyrs_cut_ok[[1]][[1]])
plotCun <- plotGrid(lyr = lyrs_cut_ok[[2]][[1]])
plotVal <- plotGrid(lyr = lyrs_cut_ok[[3]][[1]])

plotAll <- grid.arrange(plotAnt, plotCun, plotVal, ncol=3)

Aquí podemos observar la precipitación para el mes de enero en cada departamento:

Muchas gracias a Jeison Mesa por la ayuda brindada en la construcción de este tutorial.

Función Zonal Histogram de ArcGIS en R

En este tutorial explicare una aplicación de la función “zonal histogram” de ArcGIS en R, el objetivo, será calcular el área de las coberturas para los municipios del Valle del Cauca a partir de la capa de cobertura de la NASA y la división político administrativa del SIGOT.
Los datos a utilizar, en coordenadas geográficas WGS 1984, son los siguientes (disponibles aquí):
1.      Cobertura de la tierra reclasificada de la NASA.
2.      Municipios del Valle del Cauca
A continuación los pasos:
1.      Cargamos las librerías a utilizar, en caso que usted no las tenga instalada podría utilizar install.packages(“nameLibrary”)

# Cargamos librerias
require(raster)
require(rgdal)
require(foreign)
require(ggplot2)
require(dplyr)
2.      Cargamos los archivos, que son el raster, la tabla de la leyenda del raster, y el shapefile de municipios.

path <- "E:/R/_tutoriales/_tabulateArea/_blog/_data"
cob <- raster(paste0(path, "/_raster/cov_eluglcde.tif"))
adm1 <- shapefile(paste0(path, "/_shp/mpiosValle.shp"))
leg <- read.dbf(paste0(path, "/_raster/cov_eluglcde.tif.vat.dbf"))

3.      Ploteamos el shapefile para observar los límites de los municipios.

ggplot() + geom_polygon(data=adm1,  aes(x=long, y=lat, group=group), fill="cadetblue", color="grey") + coord_equal()
4.      Creamos la función de zonal histogram, pues esta función aún no existe como tal en alguna librería.

tabFunc                            <- function(indx, extracted, region, regname) {
dat <- as.data.frame(table(extracted[[indx]]))
dat$name <- region[[regname]][[indx]]
return(dat)
}
5.      Ejecutamos la función creada anteriormente

ext <- extract(cob, adm1, method='simple')#extract valores
tabs <- lapply(seq(ext), tabFunc, ext, adm1, "NOM_MUNICI")
df <- do.call(rbind, tabs)#unimos todas las tablas en una sola
6.      Calculamos el porcentaje de cobertura por municipio y las hectáreas, teniendo en cuenta que cada píxel tiene un área aproximada de 0.20833 hectáreas

df  %
group_by(name)%>%
mutate(Porcentaje = (Freq/sum(Freq))*100)%>%
ungroup()%>%
mutate(HA = Freq * 0.2083333)

7.      Editamos los encabezados y unimos con la tabla de leyenda para conocer cada valor de celda a que categoría corresponde ej. El valor 9 corresponde a áreas urbanas.

colnames(df)      <- c("Category", "Count", "Municipio", "Porcentaje", "Has")
df_all <- merge(df, leg, by.x = "Category", by.y = "Value")
df_all <- df_all[,c("Category", "Count.x", "Municipio", "Porcentaje", "Has", "ELU_GLC_De")]
df_all <- df_all[order(df_all$Municipio, decreasing = FALSE), ]
colnames(df_all) <- c("Category", "Count", "Municipio", "Porcentaje", "Has", "Count2", "Cobertura")
      8.       Creamos gráfica para tres municipios del Valle del Cauca, esta gráfica representa el porcentaje de cada cobertura en el municipio según el raster de cobertura de la NASA para el año 2015.

mpios <- filter(df_all, Municipio %in% c("YUMBO", "CALI", "PALMIRA"))
gg <- ggplot(mpios, aes(x = factor(Municipio), y = Porcentaje, group = Cobertura)) +
geom_bar(aes(fill = Cobertura), stat = "identity", position = "dodge") +
theme(legend.position = "bottom") + xlab("Municipio")
ggsave(plot = gg, paste0(path, "/cob_mpios_porc.png"), width = 9, height = 8, units = "in")
gg
dev.off()


9.      Guardamos la tabla como un archivo csv

write.csv(df_all, paste0(path, "/_table/zonalHistogram.csv"))

Gracias a Jeison Mesa por su colaboración en la realización de este tutorial!