Relieve
shp_mpio_aru = st_read('../../arauca/81_ARAUCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp')
## Reading layer `MGN_MPIO_POLITICO' from data source
## `C:\Users\ADMIN\Documentos\CARLOS\UNAL\Asignaturas\Geomatica\arauca\81_ARAUCA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 10 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -72.36662 ymin: 6.036228 xmax: -69.42756 ymax: 7.104381
## Geodetic CRS: WGS 84
elev = raster("../Clase 5 - DEM - 181121/arauca_dem_z10.tif")
DEM
- Recorter DEM a la geometria del departamento
elev_crop = crop(elev, shp_mpio_aru)
elev_crop = mask(elev_crop, shp_mpio_aru)
- Pendientes, Aspecto, y Colinas
slope = terrain(elev_crop, opt='slope', unit='degrees')
aspect = terrain(elev_crop, opt='aspect', unit='degrees')
hill = hillShade(slope = slope, aspect = aspect, angle = 40, direction = 270)
- convirtiendo raster a tablas
elev_crop_df = as.data.frame(elev_crop, xy = TRUE)
slope_df = as.data.frame(slope, xy = TRUE)
aspect_df = as.data.frame(aspect, xy = TRUE)
hill_df = as.data.frame(hill, xy = TRUE)
quantile(slope_df$slope, na.rm=T)
## 0% 25% 50% 75% 100%
## 0.0000000 0.4018827 0.7211327 1.4974133 82.4996909
ggplot()+
geom_tile(data = elev_crop_df, aes(x,y, fill = arauca_dem_z10), alpha = 0.75)+
scale_fill_gradientn(colors = terrain.colors(15), na.value = "transparent")+
labs(x = 'Longitud', y = 'Latitud', fill = 'Altitud [msnm]')+
ggnewscale::new_scale_fill()+
geom_tile(data = aspect_df, aes(x,y, fill = aspect), alpha = 0.5)+
scale_fill_gradientn(colors = brewer.pal(n=8, name = 'Greys'),
na.value="transparent",
name = 'Altitud [msnm]', guide = 'none')+
geom_sf(data = shp_mpio_aru, fill = NA, size = 1)+
geom_sf_text(data = shp_mpio_aru, aes(label = MPIO_CNMBR), size = 3)+
theme_bw()+
theme(legend.position = 'bottom',
legend.key.width = unit(3, 'line'),
text = element_text(size = 15))
# elev_crop = crop(elev, mpio_sp)
# elev_crop = mask(elev_crop, mpio_sp)
# plot(elev_crop, main="Elevación recortada")
# plot(mpio_sp, add=TRUE)
# text(coordinates(mpio_sp), col="black", cex=0.6,
# labels=as.character(mpio$MPIO_CNMBR))
Precipitación
Descargar y descomprimir los archivos *.tif de CHIRPS-2.0 mensual de los años comprendidos entre 2010-2019
# # DATOS CHIRPS ANUALES 2010-2019
#
# url = 'https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_monthly/tifs/'
# meses = c(paste0(rep(0,9), 1:9), 10:12)
# path_ch = 'chirps/meses/'
#
# # todos los meses en cada año
#
# library(parallel)
#
# nc = detectCores()
# cl = makeCluster(nc)
#
# clusterExport(cl, c('url'))
# clusterEvalQ(cl, {ls()})
#
# for(year in 2011:2020){
# print(year)
# path_year = paste0(path_ch, year, '/')
# if(!dir.exists(path_year)){dir.create(path_year)}
#
# parLapply(cl, meses, function(mes, year, path_y){
# name_tif = paste0('chirps-v2.0.', year, '.', mes, '.tif.gz')
# url_tif = paste0(url, name_tif)
#
# download.file(url_tif, paste0(path_y, name_tif))
# R.utils::gunzip(paste0(path_y, name_tif), remove = TRUE, overwrite = TRUE)
# }, year = year, path_y = path_year)
# }
#
# stopCluster(cl)
Recortar cada archivo *.tif de precipitación con el polygono del departamento y convertirlo en una tabla
#
# years = 2011:2020
# path_meses = 'chirps/meses/'
# files_tif = dir(path_meses, pattern = '.tif', recursive = T)
# meses = c(paste0(rep(0,9), 1:9), 10:12)
#
# cap_chirps = list()
#
# nc = detectCores()
# cl = makeCluster(nc-1)
# clusterExport(cl, c('path_meses', 'shp_mpio_aru'))
# clusterEvalQ(cl, {
# library(sf)
# library(raster)
# library(stringr)
# })
#
# cap_chirps = parLapply(cl, files_tif, function(mes_tif){
# chirps_tot = raster(paste0(path_meses, mes_tif))
# chirps_crop = crop(chirps_tot, shp_mpio_aru)
# chirps_crop = mask(chirps_crop, shp_mpio_aru)
# chirps_crop
# })
# stopCluster(cl)
# names(cap_chirps) = paste0('m',str_sub(files_tif, 18, -5))
# cap_chirps_join = stack(cap_chirps)
# cap_chirps_df = cap_chirps_join %>%
# as.data.frame(xy=TRUE)
# cap_chirps_df_l = cap_chirps_df %>%
# melt(c('x','y')) %>%
# mutate(year = rep(2011:2020, each = nrow(cap_chirps_df)*12),
# mes = rep(1:12, each = nrow(cap_chirps_df), 10))
cap_chirps_df_l = readRDS('cap_chirps_df_l.rds')
cap_chirps_df_l
Precipitación total por año, mes y mes-año
pcp_tot_year_mes = cap_chirps_df_l %>%
group_by(year = factor(year), mes = factor(mes)) %>%
summarise(mpcp = mean(value, na.rm=T))
pcp_tot_mes = cap_chirps_df_l %>%
group_by(x,y,mes) %>%
summarise(spcp = sum(value))
pcp_tot_year = cap_chirps_df_l %>%
group_by(x,y,year) %>%
summarise(spcp = sum(value))
pcp_tot_year$gpcp = cut(pcp_tot_year$spcp/100, 9)
Distribución espacial de la precipitación total anual en el departamento de arauca
ggplot()+
geom_raster(data = pcp_tot_year,
aes(x,y, fill = spcp))+
geom_sf(data = shp_mpio_aru, fill=NA)+
geom_sf_text(data = shp_mpio_aru, aes(label = MPIO_CNMBR), size = 3)+
facet_wrap(~year, ncol=2)+
scale_fill_distiller(palette = "Blues", direction = 1, na.value="transparent")+
# scale_fill_continuous(na.value="transparent")+
labs(x='Longitud', y='Latitud', fill='Precipitación [mm]')+
theme_bw()+
theme(text = element_text(size = 15))
Distribución mensual de la precipitación en el departamento de Arauca
pcp_tot_year_mes %>%
ggplot()+
aes(mes, mpcp, fill = year)+
geom_col(position = 'dodge')+
scale_fill_brewer(palette = 'Paired')+
labs(x='Mes', y='Precipitación [mm]', fill='Año')+
theme_bw()+
theme(legend.position = 'bottom',
legend.key.width = unit(3, 'line'),
text = element_text(size = 15))
Temperatura
# tmin = climateR::getTerraClim(AOI = shp_mpio_aru, param = 'tmin',
# startDate = "2011-01-01",
# endDate = "2019-12-31")
# tmax = climateR::getTerraClim(AOI = shp_mpio_aru, param = 'tmax',
# startDate = "2011-01-01",
# endDate = "2019-12-31")
#
# tmin_df = as.data.frame(tmin$terraclim_tmin, xy=TRUE)
# tmin_crop = crop(tmin$terraclim_tmin, shp_mpio_aru)
# tmin_crop = mask(tmin_crop, shp_mpio_aru)
# tmin_crop_df = as.data.frame(tmin_crop, xy=T)
# tmin_crop_df_l = tmin_crop_df %>%
# melt(c('x','y')) %>%
# mutate(year = rep(2011:2019, each = nrow(tmin_crop_df)*12),
# mes = rep(1:12, each = nrow(tmin_crop_df),9))
#
# tmax_df = as.data.frame(tmax$terraclim_tmax, xy=TRUE)
# tmax_crop = crop(tmax$terraclim_tmax, shp_mpio_aru)
# tmax_crop = mask(tmax_crop, shp_mpio_aru)
# tmax_crop_df = as.data.frame(tmax_crop, xy=T)
# tmax_crop_df_l = tmax_crop_df %>%
# melt(c('x','y')) %>%
# mutate(year = rep(2011:2019, each = nrow(tmax_crop_df)*12),
# mes = rep(1:12, each = nrow(tmax_crop_df),9))
#
# saveRDS(list(tmin_crop_df_l, tmax_crop_df_l), 'temp.rds')
temp = readRDS('temp.rds')
temp[[1]]$tipo = 'min'
temp[[2]]$tipo = 'max'
temp_tot = rbind(temp[[1]], temp[[2]])
temp_prom_year = temp_tot %>%
dcast('x+y+variable+year+mes~tipo') %>%
group_by(x,y, year) %>%
summarise(med_tipo = mean((min+max)/2))
temp_rang_mes = temp_tot %>%
group_by(year, mes, tipo) %>%
summarise(med = mean(value, na.rm=T))
temp_rang_mes %>%
ggplot()+
aes(factor(mes), med, color = tipo)+
geom_line(aes(group = interaction(tipo, year)), size = 3)+
scale_color_manual(values = c("#FC8D59", "#91BFDB"))+
labs(x='Mes', y='Tempratura [°C]', color = NULL)+
theme_bw()+
theme(legend.position = 'bottom',
legend.key.width = unit(3, 'line'),
text = element_text(size = 15))
ggplot()+
geom_raster(data = temp_prom_year,
aes(x,y, fill = med_tipo))+
geom_sf(data = shp_mpio_aru, fill=NA, color = '#000000', size = 1)+
geom_sf_text(data = shp_mpio_aru, aes(label = MPIO_CNMBR), size = 3)+
facet_wrap(~year)+
labs(x='Longitud', y='Latitud', fill='Tempratura\npromedio [°C]')+
scale_fill_distiller(palette = "RdYlBu",
na.value="transparent", breaks = seq(0,30,5))+
theme_bw()+
theme(legend.position = 'bottom',
legend.key.width = unit(3, 'line'),
text = element_text(size = 15))
