# library(AOI)
library(climateR)
# library(climate)
library(sf)
## Warning: package 'sf' was built under R version 4.1.2
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
library(raster)
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.1.2
library(rasterVis)
## Loading required package: lattice
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.0 v forcats 0.5.1
## Warning: package 'readr' was built under R version 4.1.2
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()
library(leaflet)
library(RColorBrewer)
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
param_meta$terraclim
## common.name call description units
## 1 aet aet Actual Evapotranspiration mm
## 2 water_deficit def Climatic Water Deficit mm
## 3 palmer PDSI Palmer Drought Severity Index
## 4 pet pet Reference Evapotranspiration mm
## 5 prcp ppt Accumulated Precipitation mm
## 6 q q Runoff mm
## 7 soilm soil Soil Moisture at End of Month mm
## 8 srad srad Downward Shortwave Radiation Flux at the Surface W/m2
## 9 swe swe Snow Water Equivalent at End of Month mm
## 10 tmax tmax Maximum 2-m Temperature C
## 11 tmin tmin Minimum 2-m Temperature C
## 12 vp vap 2-m Vapor Pressure kPa
## 13 vpd vpd Mean Vapor Pressure Deficit kPa
## 14 wind ws Wind Speed at 10-m m/s
tc_prcp = getTerraClim(shp_mpio_aru, param = "prcp",
startDate = "2019-01-01",
endDate = "2019-12-31")
## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on
plot(tc_prcp$terraclim_prcp, col = brewer.pal(n = 10, name = 'RdYlBu'))
tc_tmp <- tc_prcp$terraclim_prcp$X2019.02
tc_tmp_crop = crop(tc_tmp, shp_mpio_aru)
tc_tmp_crop = mask(tc_tmp_crop, shp_mpio_aru)
# pal <- colorNumeric(c("red", "orange", "#fcc000",
# "yellow", "cyan", "blue", "#3240cd"),
# values(tc_tmp_crop),
# na.color = "transparent")
pal = colorNumeric(palette = brewer.pal(n = 7, name = 'RdYlBu'),
domain = values(tc_tmp_crop),
na.color = "transparent")
# plot(rep(1:7, 2), rep(1:2, each=7), pch = 16, cex = 5,
# col=c(
# c("red", "orange", "#fcc000","yellow", "cyan", "blue", "#3240cd"),
# brewer.pal(n = 7, name = 'RdYlBu')
# ))
leaflet(shp_mpio_aru) %>%
addTiles() %>%
addPolygons(color = 'black', weight = 1, opacity = 1, fillOpacity = 0) %>%
addRasterImage(tc_tmp_crop,
colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(tc_tmp_crop),
title = "Precipitación - Febrero 2019 [mm]")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
tc_palmer = getTerraClim(shp_mpio_aru, param = "palmer",
startDate = "2019-02-01")
## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on
tc_tmp = tc_palmer$terraclim_palmer
tc_tmp
## class : RasterStack
## dimensions : 27, 71, 1917, 1 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -72.375, -69.41667, 6, 7.125 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2019.02
## min values : -4.6
## max values : -0.7
tc_tmp = tc_tmp$X2019.02
tc_tmp_crop = crop(tc_tmp, shp_mpio_aru)
tc_tmp_crop = mask(tc_tmp_crop, shp_mpio_aru)
pal = colorNumeric(palette = brewer.pal(n = 10, name = 'RdYlGn'),
domain = values(tc_tmp_crop),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(tc_tmp_crop, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(tc_tmp_crop),
title = "PDSI-Feb.2019")
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
wat_def = getTerraClimNormals(shp_mpio_aru,
param = "water_deficit",
period = "19812010", month=2)
## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on
tc_tmp = wat_def$terraclim_19812010_water_deficit
tc_tmp
## class : RasterStack
## dimensions : 27, 71, 1917, 1 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -72.375, -69.41667, 6, 7.125 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X02
## min values : 3.8
## max values : 77.1
tc_tmp = tc_tmp$X02
tc_tmp_crop = crop(tc_tmp, shp_mpio_aru)
tc_tmp_crop = mask(tc_tmp_crop, shp_mpio_aru)
pal = colorNumeric(palette = brewer.pal(n = 10, name = 'RdYlGn'),
domain = values(tc_tmp_crop),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(tc_tmp_crop, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(tc_tmp_crop),
title = "WaterDeficit-February")
chirps = getCHIRPS(shp_mpio_aru,
startDate = "2020-11-01",
endDate = "2020-11-05" )
## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on
(capas = names(chirps))
## [1] "prcp_20201101" "prcp_20201102" "prcp_20201103" "prcp_20201104"
## [5] "prcp_20201105"
cellStats(chirps, mean)
## prcp_20201101 prcp_20201102 prcp_20201103 prcp_20201104 prcp_20201105
## 44.715217 6.273913 15.860870 13.573913 15.135507
for(cp in capas){
tmp = chirps[[cp]]
cat(cp,cellStats(tmp, range),'\n')
}
## prcp_20201101 2 126
## prcp_20201102 2 71
## prcp_20201103 2 71
## prcp_20201104 2 50
## prcp_20201105 2 140
m <- leaflet() %>%
addTiles(group = "OSM (default)")
valores <- seq(150)
pal = colorNumeric(palette = brewer.pal(n = 7, name = 'RdYlBu'),
domain = valores,
na.color = "transparent")
for(cp in capas){
tmp = chirps[[cp]]
tmp_crop = crop(tmp, shp_mpio_aru)
tmp_crop = mask(tmp_crop, shp_mpio_aru)
m <- m %>%
addRasterImage(tmp_crop, colors = pal, opacity = 0.8, group=cp)
}
m <- m %>%
# Add layers controls
addLayersControl(
baseGroups = capas,
options = layersControlOptions(collapsed = FALSE)
) %>%
# Add common legend
addLegend(pal = pal, values = valores,
title = "CHIRPS - Nov.2020")
m
par(mfrow=c(2,3))
for(cp in capas){
tmp = chirps[[cp]]
tmp_crop = crop(tmp, shp_mpio_aru)
tmp_crop = mask(tmp_crop, shp_mpio_aru)
hist(tmp_crop)
}
# url = 'https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_monthly/tifs/'
# meses = c('01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12')
#
# cap_chirps = list()
# for(mes in meses){
# name_tif = paste0('chirps-v2.0.2020.', mes)
# url_tif = paste0(url, name_tif, '.tif.gz')
# download.file(url_tif, paste0("CHIRPS/", name_tif, '.tif.gz'))
# R.utils::gunzip(paste0("CHIRPS/", name_tif, '.tif.gz'),
# remove = TRUE, overwrite = TRUE)
# chirps_tot = raster(paste0("CHIRPS/", name_tif, ".tif"))
# chirps_crop = crop(chirps_tot, shp_mpio_aru)
# chirps_crop = mask(chirps_crop, shp_mpio_aru)
#
# cap_chirps[[mes]] = chirps_crop
# }
#
# cap_chirps_join = stack(cap_chirps)
files_tif = dir('CHIRPS/', pattern = '.tif')
meses = c('01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12')
cap_chirps = list()
for(mes in files_tif){
chirps_tot = raster(paste0("CHIRPS/", mes))
chirps_crop = crop(chirps_tot, shp_mpio_aru)
chirps_crop = mask(chirps_crop, shp_mpio_aru)
nombre = stringr::str_sub(mes, 12, -5)
cap_chirps[[paste0('m',nombre)]] = chirps_crop
}
cap_chirps_join = stack(cap_chirps)
capas = names(cap_chirps_join)
for(cp in capas){
tmp = cap_chirps_join[[cp]]
cat(cp,cellStats(tmp, range),'\n')
}
## m.2020.01 1.333755 114.3775
## m.2020.02 0.8956861 72.04098
## m.2020.03 2.032927 61.38647
## m.2020.04 71.14093 187.1743
## m.2020.05 75.03194 420.4474
## m.2020.06 191.5283 486.7229
## m.2020.07 163.5916 576.6873
## m.2020.08 116.1901 536.0602
## m.2020.09 124.8234 320.1478
## m.2020.10 92.34514 330.636
## m.2020.11 125.233 377.9701
## m.2020.12 7.533748 103.2016
m <- leaflet(shp_mpio_aru) %>%
addTiles(group = "OSM (default)") %>%
addPolygons(color = 'black', weight = 1, opacity = 1, fillOpacity = 0)
valores <- seq(600)
pal = colorNumeric(palette = brewer.pal(n = 7, name = 'RdYlBu'),
domain = valores,
na.color = "transparent")
for(cp in capas){
tmp = cap_chirps_join[[cp]]
m <- m %>%
addRasterImage(tmp, colors = pal,
opacity = 0.8, group=cp)
}
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
m <- m %>%
addLayersControl(
baseGroups = capas,
options = layersControlOptions(collapsed = FALSE)
) %>%
addLegend(pal = pal, values = valores,
title = "CHIRPS_v2.0 - 2020")
m