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)


Comentarios

  1. 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 extent

    he probado ha proyectar uno sobre el otro con

    projectRaster(r2,r1)

    pero sale vacio ya que son zonas distintas

    alguna solucion un saludo

    ResponderEliminar

Publicar un comentario

Entradas populares de este blog

Extracción por mascara en R

Convertir una tabla a shape en R

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