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)

Un pensamiento en “Creando un mapa en R

  1. Unknown

    hola fabio estoy aprendiendo a usar raster, tengo varios tif de parcelas del campo distintas, hay alguna forma de unirlos?los tif los cargo como si fueran raster pero el problema que las coordenadas respetan un sistema de regilla distinto, al unirlos con stack obtengo Error in compareRaster(x) : different extenthe probado ha proyectar uno sobre el otro con projectRaster(r2,r1)pero sale vacio ya que son zonas distintasalguna solucion un saludo

    Le gusta a 1 persona

    Responder

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión /  Cambiar )

Google photo

Estás comentando usando tu cuenta de Google. Cerrar sesión /  Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión /  Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión /  Cambiar )

Conectando a %s