rm(list=ls())
library(AOI)
library(climateR)
library(sf)
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(raster)
## Loading required package: sp
library(rasterVis)
## Loading required package: terra
## terra version 1.3.4
## Loading required package: lattice
## Loading required package: latticeExtra
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     intersect, src, union
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

##3. Get climate data

santander<-st_read("C:/Users/morea/OneDrive - Universidad Nacional de Colombia/Documentos/UNAL/2021 - 1S/GeomaticaBasica/MGN2017_68_SANTANDER/68_SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source 
##   `C:\Users\morea\OneDrive - Universidad Nacional de Colombia\Documentos\UNAL\2021 - 1S\GeomaticaBasica\MGN2017_68_SANTANDER\68_SANTANDER\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 87 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
## Geodetic CRS:  WGS 84
tc_prcp = getTerraClim(santander, param = "prcp", startDate = "2019-02-01")
## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on
tc_tmp<-tc_prcp[[1]]
tc_tmp
## class      : RasterStack 
## dimensions : 60, 50, 3000, 1  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -74.54167, -72.45833, 5.666667, 8.166667  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : X2019.02 
## min values :      2.4 
## max values :    101.9
library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "#fcc000","yellow", "cyan", "blue", "#3240cd"), values(tc_tmp$X2019.02), 
                    na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(tc_tmp$X2019.02 , colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(tc_tmp$X2019.02),
    title = "Rainfall-Feb.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 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
head(param_meta$terraclim)

Palmer Drought Severity Index (PDSI)

tc_palmer = getTerraClim(santander, param = "palmer", startDate = "2019-02-01")
## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on
tc_tmp <- tc_palmer[[1]]
tc_tmp
## class      : RasterStack 
## dimensions : 60, 50, 3000, 1  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -74.54167, -72.45833, 5.666667, 8.166667  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : X2019.02 
## min values :     -4.1 
## max values :      1.3
pal <- colorNumeric(c("#fc7300","orange", "yellow","#9acd32", "green"), values(tc_tmp$X2019.02), 
                    na.color = "transparent")
leaflet() %>% addTiles() %>%
  addRasterImage(tc_tmp$X2019.02, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(tc_tmp$X2019.02),
    title = "PDSI-Feb.2019")
## 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 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

###3.2 Obtaining TerraClimate normals

Water Deficit data (Indice climatico de agua)

wat_def = getTerraClimNormals(santander, param = "water_deficit", period = "19812010", month=2)
## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on
wat_def
## $terraclim_19812010_water_deficit
## class      : RasterStack 
## dimensions : 60, 50, 3000, 1  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -74.54167, -72.45833, 5.666667, 8.166667  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :  X02 
## min values :  6.8 
## max values : 78.2
tc_tmp<-wat_def[[1]]
pal <- colorNumeric(c("green", "#9acd32","yellow", "orange", 
                    "#fc7300"), values(tc_tmp$X02), 
                    na.color = "transparent")
leaflet() %>% addTiles() %>%
  addRasterImage(tc_tmp$X02, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(tc_tmp$X02),
    title = "WaterDeficit-February")
## 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 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

###3.3 Obtaining CHIRPS