Importando los datos de precipitación del producto PISCO

###Paquete para manejo de datos raster
library(raster)
###Rasterizando el archivo netcdf, por capas (Rasterlayer)
nc <-brick("pp.nc")
nc
## class      : RasterBrick 
## dimensions : 198, 133, 26334, 432  (nrow, ncol, ncell, nlayers)
## resolution : 0.09999998, 0.1  (x, y)
## extent     : -81.3, -68, -18.8, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : pp.nc 
## names      : X252.5, X253.5, X254.5, X255.5, X256.5, X257.5, X258.5, X259.5, X260.5, X261.5, X262.5, X263.5, X264.5, X265.5, X266.5, ... 
## T (months since 1960-01-01): 252.5, 683.5 (min, max)
## varname    : Prec

Cálculo del índice SPI (Standarized Precipitation Index)

Para el cálculo del índice SPI se hará usando el paquete SPI usado en múltiples estudios

SPI
## class      : RasterBrick 
## dimensions : 198, 133, 26334, 432  (nrow, ncol, ncell, nlayers)
## resolution : 0.09999998, 0.1  (x, y)
## extent     : -81.3, -68, -18.8, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : r_tmp_2023-05-30_114637_21268_82095.grd 
## names      :   layer.1,   layer.2,   layer.3,   layer.4,   layer.5,   layer.6,   layer.7,   layer.8,   layer.9,  layer.10,  layer.11,  layer.12,  layer.13,  layer.14,  layer.15, ... 
## min values :        NA,        NA, -2.909990, -3.062407, -3.722790, -4.402364, -4.642695, -2.769820, -3.562222, -2.752526, -3.702912, -3.319882, -3.175056, -2.541338, -2.921212, ... 
## max values :        NA,        NA, 2.4496528, 2.6347866, 2.4210070, 3.2746037, 2.1153402, 4.4520390, 3.3786482, 3.8333384, 3.2485798, 2.9102327, 2.7357671, 3.6833464, 3.3448609, ...
###Abriendo Shapefile de todo sudamerica
shp <- shapefile("Sudamérica.shp")

###Filtrando Peru
shp <- shp[shp$PAÍS == "Perú" ,]
shp
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -81.35515, -68.6739, -18.34855, -0.036875  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 1
## names       : PAÍS 
## value       : Perú
###Estableciendo sistema de refererencia al rasterbrick de los indices
crs(SPI) <- crs(shp)
SPI
## class      : RasterBrick 
## dimensions : 198, 133, 26334, 432  (nrow, ncol, ncell, nlayers)
## resolution : 0.09999998, 0.1  (x, y)
## extent     : -81.3, -68, -18.8, 1  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : r_tmp_2023-05-30_114637_21268_82095.grd 
## names      :   layer.1,   layer.2,   layer.3,   layer.4,   layer.5,   layer.6,   layer.7,   layer.8,   layer.9,  layer.10,  layer.11,  layer.12,  layer.13,  layer.14,  layer.15, ... 
## min values :        NA,        NA, -2.909990, -3.062407, -3.722790, -4.402364, -4.642695, -2.769820, -3.562222, -2.752526, -3.702912, -3.319882, -3.175056, -2.541338, -2.921212, ... 
## max values :        NA,        NA, 2.4496528, 2.6347866, 2.4210070, 3.2746037, 2.1153402, 4.4520390, 3.3786482, 3.8333384, 3.2485798, 2.9102327, 2.7357671, 3.6833464, 3.3448609, ...

Graficado basico de un tiempo

###Paquete para el grafciado de datos raster o netcdf
library(rasterVis)

###Visualizando los maximos o minimos del indice desde 1981 al 2016
max(maxValue(SPI),na.rm=T)
## [1] 5.934945
min(minValue(SPI),na.rm=T)
## [1] -9.282772
###Graficando una fecha cualquiera
levelplot(SPI[[3]],margin=F)

Recortando raster para el Peru

library(tidyverse)
###Recortando con el shapefile de Perú
recor <- raster::crop(SPI,shp) %>%
  raster::mask(shp)
levelplot(recor[[3]],margin=F)

Raster a array

Convirtiendo raster del SPI para Perú a arreglo

SPI_arr <- raster::as.array(recor)

str(SPI_arr)
##  num [1:183, 1:126, 1:432] NA NA NA NA NA NA NA NA NA NA ...

Extrayendo fechas del raster anterior

dates <- getZ(nc)
dates <-as.Date("1960-01-01") %m+% months(dates-0.5)
head(dates)
## [1] "1981-01-01" "1981-02-01" "1981-03-01" "1981-04-01" "1981-05-01"
## [6] "1981-06-01"

Extracción de latitud y longitud

###Extrayendo las coordenadas de latitud y longitud del raster anterior
lon <- raster::xFromCol(recor, 1:ncol(recor))
lat <- raster::yFromRow(recor, 1:nrow(recor))

###Vemos las longitudes de cada eje
c(length(lon),length(lat))
## [1] 126 183

Array invertido, longitudes deben estar primero

str(SPI_arr)
##  num [1:183, 1:126, 1:432] NA NA NA NA NA NA NA NA NA NA ...

Corrección del array

SPI_arr <- aperm(SPI_arr, c(2,1,3))
###Corregido
str(SPI_arr)
##  num [1:126, 1:183, 1:432] NA NA NA NA NA NA NA NA NA NA ...

Latitudes invertidas

head(lat)
## [1] -0.05000001 -0.15000001 -0.25000001 -0.35000001 -0.45000001 -0.55000001

Correción de la latitud

###La latitud no va de menor a mayor
lat <- rev(lat)  

#Inversión también de latitudes en el arreglo
SPI_arr <- SPI_arr[,ncol(SPI_arr):1 ,]

head(lat)
## [1] -18.25 -18.15 -18.05 -17.95 -17.85 -17.75

Correlación del SPI con índices

Descarga del índice ONI

library(rsoi)

oni <- download_oni()
oni
## # A tibble: 880 × 7
##     Year Month Date       dSST3.4    ONI ONI_month_window phase             
##    <int> <ord> <date>       <dbl>  <dbl> <chr>            <fct>             
##  1  1950 Ene.  1950-01-01   -1.62 NA     <NA>             <NA>              
##  2  1950 Feb.  1950-02-01   -1.32 -1.34  EFM              Cool Phase/La Nina
##  3  1950 Mar.  1950-03-01   -1.07 -1.17  FMA              Cool Phase/La Nina
##  4  1950 Abr.  1950-04-01   -1.11 -1.18  MAM              Cool Phase/La Nina
##  5  1950 May.  1950-05-01   -1.37 -1.07  AMJ              Cool Phase/La Nina
##  6  1950 Jun.  1950-06-01   -0.74 -0.85  MJJ              Cool Phase/La Nina
##  7  1950 Jul.  1950-07-01   -0.44 -0.533 JJA              Cool Phase/La Nina
##  8  1950 Ago.  1950-08-01   -0.42 -0.423 JAS              Neutral Phase     
##  9  1950 Set.  1950-09-01   -0.41 -0.383 ASO              Neutral Phase     
## 10  1950 Oct.  1950-10-01   -0.32 -0.443 SON              Neutral Phase     
## # ℹ 870 more rows

Filtrando por fechas a usar del índice

###Filtrando 1991-2020
oni <- oni[oni$Date>=as.Date("1981-01-01") & oni$Date<=as.Date("2016-12-31"),]
head(oni$Date)
## [1] "1981-01-01" "1981-02-01" "1981-03-01" "1981-04-01" "1981-05-01"
## [6] "1981-06-01"

Correlacionando

Matrices vacías para almacenar coeficientes de correlación y pvalores

c.matrix <- matrix(NA,length(lon),length(lat))
t.matrix <- matrix(NA,length(lon),length(lat))

Bucle para obtención de coeficiente de correlación y pvalor

Creando dataframe con coordenas , coeficientes y pvalores

grid <- expand.grid(x=lon, y=lat)
grid$corr <- as.vector(c.matrix)
grid$pval <- as.vector(t.matrix)


sig <- subset(grid[, c(1, 2, 4)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)
sig
## class       : SpatialPointsDataFrame 
## features    : 3289 
## extent      : -81.15, -68.75, -18.25, -0.25  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 3
## names       :                 x,                  y,                 pval 
## min values  : -81.1500000231194,  -18.2500000003026, 1.78638801297263e-12 
## max values  : -68.7500028899222, -0.250000011194779,  0.00999805909927293

Graficando con rasterVis

library(maps)
## 
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
## 
##     map
library(maptools)
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired during 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
library(rgdal)
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-5, (SVN revision 1199)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/ASUS/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/ASUS/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(paletteer)
countries=maps::map('world',plot=F)
countries=map2SpatialLines(countries)

mycolorkey <- list(labels=list(tri.lower = TRUE, tri.upper = TRUE,labels=c(seq(-1,1,0.2)), at=seq(-1,1,0.2)))

#png("cor2.png", width = 9, height = 10, units = 'in', res = 700)
levelplot(corr~x*y,data=grid,col.regions = rev(paletteer_c("ggthemes::Classic Red-White-Black", length(seq(-1,1,0.1))-1)),at=seq(-1,1,0.1),xlab=("Longitud (°)"),ylab=("Latitud (°)"),pretty=T,maxpixels=1e50,main="Correlación de Spearman SPI e Índice ONI \n (1981-2016) pval<0.01",colorkey=mycolorkey)  + 
  latticeExtra::layer(sp.lines(countries))+
  latticeExtra::layer(sp.points(sig, pch = 1.5, cex = 0.2, col = "black"))

#dev.off()

Graficando con ggplot

library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:raster':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(sf)
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(metR)
## 
## Attaching package: 'metR'
## The following object is masked from 'package:purrr':
## 
##     cross
peru <- st_read("Sudamérica.shp")
## Reading layer `Sudamérica' from data source 
##   `C:\Users\ASUS\Desktop\TF_Tecnicas\Sudamérica.shp' using driver `ESRI Shapefile'
## Simple feature collection with 15 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -109.4461 ymin: -58.49861 xmax: -26.24139 ymax: 12.59028
## Geodetic CRS:  WGS 84
peru <- peru[peru$PAÍS == "Perú",]
sig <- subset(grid[, c(1, 2, 4)], pval < 0.01) 

ggplot(data=peru)+
  #geom_tile(data=grid,aes(x,y,fill=corr))+
  geom_contour_fill(data=grid,aes(x,y,z=corr),breaks=seq(-1,1,0.1))+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5,face = "bold"),legend.key.height = unit(4.5,"cm"))+
  scale_fill_gradientn(colours=rev(paletteer_c("ggthemes::Red-Blue-White Diverging",length(seq(-1,1,0.1))-1)),breaks=seq(-1,1,0.1),limits=c(-1,1))+
  geom_sf(fill="transparent",col="black",lwd=0.6)+
  geom_point(data=sig,aes(x,y),size=0.00001)+
  ggtitle("Correlación de Spearman SPI-3 e Índice ONI pval<0.01")+
  xlab("Longitud")+
  ylab("Latitud")+
  guides(fill=guide_colorsteps(title="Rho"))+
  coord_sf(xlim=c(-81,-69),ylim=c(-17.55,-0.8))

Graficando con plotly

ggplotly(ggplot(data=peru)+
  geom_point(data=sig,aes(x,y),size=0.015)+
  geom_tile(data=grid,aes(x,y,fill=corr),breaks=seq(-1,1,0.1))+theme_bw()+
  theme(plot.title = element_text(hjust = 0.5,face = "bold"))+
  scale_fill_gradientn(colours=rev(paletteer_c("ggthemes::Red-Blue-White Diverging", 30)),breaks=seq(-1,1,0.2),limits=c(-1,1))+
  geom_sf(fill="transparent",col="black",lwd=0.5)+
  ggtitle("Correlación de Spearman SPI-3 e Índice ONI \n (1983-2016) pval<0.01")+
  xlab("Longitud")+
  ylab("Latitud")+
  coord_sf(xlim=c(-81,-69),ylim=c(-17.55,-0.5)))

Bucle para automatizar las gráficas de otros índices

Establecer nombres a la lista

names(graph) <- name_index
names(graph)
## [1] "ONI"  "SOI"  "AAO"  "AO"   "PDO"  "DMI"  "MEI"  "NPGO"

Gráficas

graph$ONI

graph$SOI

graph$AAO

graph$AO

graph$PDO

graph$DMI

graph$MEI

graph$NPGO