Shapefiles (polígonos, puntos, lineas)

Cargamos paquetes a utilizar

El paquete más utilizado para trabajar con shapefiles es el sf, entonces lo primero que haremos será cargar ese paquete

library(sf)

Cargando capas

La función para cargar capas de tipo shapefile es read_sf cada vez que se use una función de Rstudio se debe poner Insert (botón en la parte superior derecha) Para mostrar como tabla el shape que se lee se tiene que asignar un nombre. Se pone “nombre que desees” <- read_sf (el archivo shapefile)

Departamentos <- read_sf("Basica/limites_politicos/BAS_LIM_DEPARTAMENTO.shp")

Graficando mapas

usando sf

Para graficar se usa la funcion plot. Cuando solo se utiliza “plot(Nombre del archivo shipe)” se graficara cada columna de la tabla. Ahora para escoger que se desea graficar se usa corchetes. “plot(nombre del archivo shipe[nombre de la columna a graficar])”

Para saber el nombre de las columnas se usa la función colnames. “colnames(nombre del archivo shipe)”

colnames(Departamentos)
## [1] "NOMBDEP"    "COUNT"      "FIRST_IDDP" "HECTARES"   "geometry"

Entonces voy a graficar la columna HECTARES El grafico aparecera con su leyenda. Para la ubicación de la leyenda se coloca la funcion key.pos. Esta puede aparece: Arriba = 4 Abajo = 1 Derecha = 3 Izquierda = 2

plot(Departamentos["HECTARES"], key.pos = 1)

Para escoger mas de una columna a graficar. Se usa c cuando se usa mas de dos elementos. “plot(archivo shipe[c(“nombre columna 1”,“nombre columna 2”, ….)], key.pos = 1)”

plot(Departamentos[c("HECTARES","COUNT")], key.pos = 1)

Graficar solo geometría

Si se desea graficar solo la geometria del doc. usar la funcion st_geometry. “plot(st_geometry(nombre del archivo))”

plot(st_geometry(Departamentos))

Color y ejes de gráficos

Para poner las caracteristicas de color y ejes. Usar la siguiente función. “plot(st_geometry(nombre del archivo), col = sf.colors(cuantos colores deseas, categorical = TRUE), border = ‘nombre del color de borde’, axes = TRUE)

NOTA: la funcion axes es para especificar si se desea poner ejes. (TRUE O FALSE)

plot(st_geometry(Departamentos), col = sf.colors(24, categorical = TRUE), border = 'black', 
axes = TRUE)

Controlar el número de gráficos

Usar la función max.plot. “plot(archivo shipe, max.plot = # de graficos que se desea”

plot(Departamentos, max.plot=2)

Para cambiar de unidades en este caso de hectarias a km2

Primero se carga dyplyr. “library(dyplyr)” Usar mutate. “shipe %>% mutate(unidad a la que se quiere cambiar = unidad a cambiar/equivalente)”

library(dplyr)
Departamentos <- Departamentos %>% mutate(km2 = HECTARES/100)

Para controlar los rangos de los ejes

Se usa las funciones de breaks y nbreaks. “plot(shipe[“columna”], breaks = c(intervalos))”

plot(Departamentos["HECTARES"], breaks = c(0,1000000,3000000,5000000,7000000,9000000))

Mapas interactivos en pagina html

Cargar el paquete leaflet, asegurarse que la projección del mapa este en latitud-longitud, en este caso al ver el st_crs, vemos que esta en metros:

st_crs(Departamentos)
## Coordinate Reference System:
##   EPSG: 32718 
##   proj4string: "+proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs"

Para transformar en latitud y longitud usar la funcion st_transform

library(leaflet)

Departamentos <- Departamentos %>% st_transform (crs = 4326)
leaflet(Departamentos) %>% addTiles() %>% addPolygons()

Si solo selecciono algunos departamentos

DepartamentosSeleccionados <- Departamentos %>% dplyr::filter(NOMBDEP %in% c("AMAZONAS", "HUANUCO"))

leaflet(DepartamentosSeleccionados) %>% addTiles() %>% addPolygons()

Seleccionar aquellos departamento con mas de 9000000 hectareas. Si se desea agregar los nombres usar “addPolygons(popup = ~columna del nombre)”

DepartamentosGrandes <- Departamentos %>% dplyr::filter(HECTARES >= 9000000)

leaflet(DepartamentosGrandes) %>% addTiles() %>% addPolygons(popup = ~NOMBDEP)

Recorte de mapas

Rasters (grillas)

Primero cargo los paquetes

library(raster)
library(rgdal)

Algunas librerias extra

  • Para visualizar el paquete rasterVis

Para bajar datos con el paquete raster

Se usa función getData

a continución bajaremos una capa de Precipitación

Prec <- getData("worldclim", var = "prec", res = 10)

Para graficar con rasterVis

library(rasterVis)
plot(Prec)

Para cortar uso crop

Prec_Peru <- crop(Prec, Departamentos)

y ahora lo grafico

plot(Prec_Peru)

Para un corte más fino uso la función mask

Prec_Peru2 <- mask(Prec_Peru, Departamentos)
plot(Prec_Peru2)

levelplot(Prec_Peru2)

Si yo quiero subsettear un stack puedo usar dos parentesis cuadrados

quiero ver la precipitación de julio

Julio <- Prec_Peru2[[7]]
plot(Julio)

Invierno <- Prec_Peru2[[c(7, 8, 9)]]
levelplot(Invierno)

Supongamos que quiero sumar todas mis capas y generar una capa de pp anual

PPAnual<- do.call("sum", unstack(Prec_Peru2))
plot(PPAnual)

la precipitacion promedio de todo perú es

cellStats(PPAnual, "mean")
## [1] 1564.149

Estudiemos el cambio climático

Haremos lo mismo pero al futuro

Futuro_CC_2070_85 <- getData("CMIP5", var = "prec", res = 10, year= 70, model= "CC", rcp = 85)
plot(Futuro_CC_2070_85)

Enero <- Futuro_CC_2070_85[[1]]
plot(Enero)

Futuro_PPeru <- crop(Futuro_CC_2070_85, Departamentos)
plot(Futuro_PPeru)

Futuro_PPeru2 <- mask(Futuro_PPeru, Departamentos)
plot(Futuro_PPeru2)

levelplot(Futuro_PPeru2)

PPFuturo<- do.call("sum", unstack(Futuro_PPeru2))
plot(PPFuturo)

Diferencia <- PPFuturo - PPAnual
plot(Diferencia)