rm(list=ls())
library(AOI)
library(climateR)
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(raster)
## Loading required package: sp
library(rasterVis)
## Loading required package: lattice
library(dplyr)
## 
## Attaching package: 'dplyr'
## 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
library(sf)
(narino <- sf::st_read("/cloud/project/MunicipiosNarino.shp"))
## Reading layer `MunicipiosNarino' from data source 
##   `/cloud/project/MunicipiosNarino.shp' using driver `ESRI Shapefile'
## Simple feature collection with 66 features and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS:  WGS 84
## Simple feature collection with 66 features and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS:  WGS 84
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR
## 1          52        083      BELÉN
## 2          52        110    BUESACO
## 3          52        203      COLÓN
## 4          52        480     NARIÑO
## 5          52        506     OSPINA
## 6          52        720    SAPUYES
## 7          52        786  TAMINANGO
## 8          52        788     TANGUA
## 9          52        240  CHACHAGÜÍ
## 10         52        254   EL PEÑOL
##                                                     MPIO_CRSLC MPIO_NAREA
## 1                            Ordenanza 53 Noviembre 29 de 1985   41.84541
## 2                                                         1899  635.96083
## 3                                         Ordenanza 37 de 1921   61.75053
## 4  Ordenanza 027 de Noviembre. 29 de 1999. Decreto 0312 del 24   25.31281
## 5                                         Ordenanza 50 de 1865   64.84321
## 6                                                         1849  115.54851
## 7                                                         1834  234.65783
## 8                                        Ordenanza 103 de 1874  217.95977
## 9                            Ordenanza 20 Noviembre 24 de 1992  146.27176
## 10                        Ordenanza 036 de Diciembre 7 de 1998  119.85744
##    MPIO_CCNCT MPIO_NANO DPTO_CNMBR  SHAPE_AREA SHAPE_LEN ORIG_FID
## 1       52083      2020     NARIÑO 0.003391678 0.3732840        0
## 2       52110      2020     NARIÑO 0.051533090 1.2292312        1
## 3       52203      2020     NARIÑO 0.005005108 0.4592866        2
## 4       52480      2020     NARIÑO 0.002050175 0.2642048        3
## 5       52506      2020     NARIÑO 0.005249269 0.3371496        4
## 6       52720      2020     NARIÑO 0.009351438 0.6599792        5
## 7       52786      2020     NARIÑO 0.019009395 0.6601636        6
## 8       52788      2020     NARIÑO 0.017652117 0.7808421        7
## 9       52240      2020     NARIÑO 0.011849554 0.7373282        8
## 10      52254      2020     NARIÑO 0.009707107 0.5550189        9
##                          geometry
## 1  POLYGON ((-77.07227 1.63422...
## 2  POLYGON ((-77.23516 1.45240...
## 3  POLYGON ((-77.04473 1.67173...
## 4  POLYGON ((-77.34282 1.31465...
## 5  POLYGON ((-77.55776 1.07006...
## 6  POLYGON ((-77.71499 1.0915,...
## 7  POLYGON ((-77.32644 1.67981...
## 8  POLYGON ((-77.36152 1.19568...
## 9  POLYGON ((-77.30295 1.51777...
## 10 POLYGON ((-77.39239 1.60127...
narino %>% dplyr::select( MPIO_CCNCT, MPIO_CNMBR) -> narino
narino
## Simple feature collection with 66 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS:  WGS 84
## First 10 features:
##    MPIO_CCNCT MPIO_CNMBR                       geometry
## 1       52083      BELÉN POLYGON ((-77.07227 1.63422...
## 2       52110    BUESACO POLYGON ((-77.23516 1.45240...
## 3       52203      COLÓN POLYGON ((-77.04473 1.67173...
## 4       52480     NARIÑO POLYGON ((-77.34282 1.31465...
## 5       52506     OSPINA POLYGON ((-77.55776 1.07006...
## 6       52720    SAPUYES POLYGON ((-77.71499 1.0915,...
## 7       52786  TAMINANGO POLYGON ((-77.32644 1.67981...
## 8       52788     TANGUA POLYGON ((-77.36152 1.19568...
## 9       52240  CHACHAGÜÍ POLYGON ((-77.30295 1.51777...
## 10      52254   EL PEÑOL POLYGON ((-77.39239 1.60127...
rename(narino, MPIO_CCDGO = MPIO_CCNCT)
## Simple feature collection with 66 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS:  WGS 84
## First 10 features:
##    MPIO_CCDGO MPIO_CNMBR                       geometry
## 1       52083      BELÉN POLYGON ((-77.07227 1.63422...
## 2       52110    BUESACO POLYGON ((-77.23516 1.45240...
## 3       52203      COLÓN POLYGON ((-77.04473 1.67173...
## 4       52480     NARIÑO POLYGON ((-77.34282 1.31465...
## 5       52506     OSPINA POLYGON ((-77.55776 1.07006...
## 6       52720    SAPUYES POLYGON ((-77.71499 1.0915,...
## 7       52786  TAMINANGO POLYGON ((-77.32644 1.67981...
## 8       52788     TANGUA POLYGON ((-77.36152 1.19568...
## 9       52240  CHACHAGÜÍ POLYGON ((-77.30295 1.51777...
## 10      52254   EL PEÑOL POLYGON ((-77.39239 1.60127...
library(climateR)
tc_prcp = getTerraClim(narino, 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 : 57, 53, 3021, 1  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -79.04167, -76.83333, 0.3333333, 2.708333  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : X2019.02 
## min values :     61.8 
## max values :    580.5
library(leaflet)
library(RColorBrewer)
library(raster)
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 colors(.): Some values were outside the color scale and will be
## treated as NA
(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
library(climateR)
library(s2)
tc_prcp = getTerraClim(narino, 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 : 57, 53, 3021, 1  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -79.04167, -76.83333, 0.3333333, 2.708333  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : X2019.02 
## min values :     61.8 
## max values :    580.5
library(leaflet) 
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 colors(.): Some values were outside the color scale and will be
## treated as NA
library(climateR)
wat_def = getTerraClimNormals(narino, 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 : 57, 53, 3021, 1  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -79.04167, -76.83333, 0.3333333, 2.708333  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :  X02 
## min values :    0 
## max values : 31.8
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 colors(.): Some values were outside the color scale and will be
## treated as NA

3.2 Obtaining climatic variables

Precipitacion

oneyear = seq(1,12)
library(climateR)
library(AOI)
precip = getTerraClimNormals(narino, param = "prcp", period = "19812010", month=oneyear)
## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on
precip
## $terraclim_19812010_prcp
## class      : RasterStack 
## dimensions : 57, 53, 3021, 12  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -79.04167, -76.83333, 0.3333333, 2.708333  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :   X01,   X02,   X03,   X04,   X05,   X06,   X07,   X08,   X09,   X10,   X11,   X12 
## min values :  40.0,  53.5,  66.7,  89.4,  44.5,  23.3,  12.3,  12.8,  27.4,  59.0,  68.6,  53.5 
## max values : 443.8, 404.0, 408.1, 533.4, 691.7, 548.0, 478.5, 439.5, 552.8, 559.4, 513.1, 522.9
capas <-names(precip)
par(mfcol=c(12,12), mfrow=c(4,3), oma=c(1,1,0,0),
    mar=c(1,1,1,0), tcl=-0.1, mgp=c(0,0,0))
for (penta in capas) {
  hist(precip[[penta]], ylab=NA, cex.axis=0.5, font.main=1, cex.main=0.8)
}

year_precip <- sum(precip$terraclim_19812010_prcp)
year_precip
## class      : RasterLayer 
## dimensions : 57, 53, 3021  (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -79.04167, -76.83333, 0.3333333, 2.708333  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 608.7, 5597.1  (min, max)
ppal <- colorNumeric(c("orange", "yellow", "#a1d99b", "#deebf7", "#9ecae1","#3182bd", "darkblue"), values(year_precip$layer), 
                    na.color = "transparent")
leaflet() %>% addTiles() %>%
  addRasterImage(year_precip, colors = ppal, opacity = 0.8) %>%
  addLegend(pal = ppal, values = values(year_precip$layer),
    title = "Annual Precipitation")
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
writeRaster(tc_tmp, filename="precip", format="GTiff", overwrite=TRUE)

Bibliografia Lizarazo, I., 2021. A tutorial on getting global TerraClimate and CHIRPS data from R. Available at https://rpubs.com/ials2un/trrclmt.