Tendencias

Instalamos el paquete

install.packages("trend")

Carga de datos

Se realiza la carga de datos de precipitación de 4 estaciones:

# Llamado del paquete
library(trend)
df <- read.csv("Precipitacion.csv")
head(df)

           X CHQ HNE  HMO PZC
1 2003-01-01 0.0 0.0 20.7 9.6
2 2003-01-02 4.7 0.8  0.9 0.0
3 2003-01-03 1.2 0.0  0.0 0.0
4 2003-01-04 0.0 0.0  0.0 0.0
5 2003-01-05 0.0 0.0  0.0 0.0
6 2003-01-06 0.0 0.0  1.7 0.2

Realizar el análisis de tendencia para cada estación meteorológica:

trend1 <- mk.test(df$CHQ)
trend1

    Mann-Kendall trend test

data:  df$CHQ
z = 2.6653, n = 5114, p-value = 0.007693
alternative hypothesis: true S is not equal to 0
sample estimates:
           S         varS          tau 
2.638330e+05 9.798906e+09 2.822036e-02 

p > 0.05: Decimos que no existe tendencia en los datos (se acepta la hipótesis nula).

p <0.05: Decimos que existe tendencia en los datos (se acepta la hipótesis alterna).

Si z > 0: Decimos que la tendencia es creciente.

Si z < 0: Decimos que la tendencia es decreciente.

Por lo tanto la estación Chuquibambilla presenta tendencia (p = 0.007693) y posee una tendencia creciente (z = 2.6653).

lista1 <- lapply(df[-1], mk.test)
lista1

$CHQ

    Mann-Kendall trend test

data:  X[[i]]
z = 2.6653, n = 5114, p-value = 0.007693
alternative hypothesis: true S is not equal to 0
sample estimates:
           S         varS          tau 
2.638330e+05 9.798906e+09 2.822036e-02 


$HNE

    Mann-Kendall trend test

data:  X[[i]]
z = -0.33184, n = 5114, p-value = 0.74
alternative hypothesis: true S is not equal to 0
sample estimates:
            S          varS           tau 
-3.247100e+04  9.574432e+09 -3.523093e-03 


$HMO

    Mann-Kendall trend test

data:  X[[i]]
z = -0.21713, n = 5114, p-value = 0.8281
alternative hypothesis: true S is not equal to 0
sample estimates:
            S          varS           tau 
-2.121800e+04  9.548132e+09 -2.305656e-03 


$PZC

    Mann-Kendall trend test

data:  X[[i]]
z = 1.7087, n = 5114, p-value = 0.08751
alternative hypothesis: true S is not equal to 0
sample estimates:
           S         varS          tau 
1.527200e+05 7.988631e+09 1.843629e-02 

Pasando de datos diarios a mensuales

library(hydroTSM)
library(xts)
library(zoo)

dfm <- as.xts(
  x = read.csv.zoo(
    "Precipitacion.csv"))
dfm <- daily2monthly(dfm, FUN = sum)
dfy <- monthly2annual(dfm, FUN = sum)

Graficación

lattice::xyplot(dfm)

Análisis de tendencia espacial

Prueba de tendencia de Mann-Kendall Se trabajará con rasters de temperatura mínima (12 rasters que representan a 12 meses del año) y se creará una lista de rasters (stack) para luego obtener otro raster de tipo lista que posea los valores estadísticos de mann-kendall (tendencai y pvalor).

library(raster)

Loading required package: sp

library(Kendall)
library(terra)

terra 1.7.39

Attaching package: 'terra'

The following object is masked from 'package:hydroTSM':

    extract

The following object is masked from 'package:zoo':

    time<-

source("Piuray/MKraster.R")

Cargando los archivos raster

nombres <- list.files(full.names = TRUE,
                      pattern = ".tif")
nombres

# Creando el raster combinado:
combinacion <- stack(nombres)
# Guardando el raster combinado
writeRaster(combinacion, "stacked.tif")

Se recomienda limpiar la memoria luego de este paso a fin de evitar saturación de datos (pero volver a cargar el archivo mk raster pues contiene el algoritmo facilitador).

# CArga del raster guardado:
img <- stack("Piuray/stacked.tif")
plot(img)

Cálculo de tendencia y significación estadística

mk <- MKraster(img, "both")

[1] "Start MKraster: 2023-07-16 02:19:55.711421"
[1] "Loading parameters"
[1] "Done loading parameters"
[1] "Initiating loop operation"
[1] "Populating raster brick"
[1] "Ending MKraster on 2023-07-16 02:20:37.15411"

Visualización de resultados

plot(mk)

Se aprecia que la tendencia es negativa (valores trend inferiores a 0) y se ve que los pvalores no se muestran significativos.

library(gridExtra)
library(rasterVis)

Loading required package: lattice

p1 <- levelplot(
  mk[[1]],names.attr = "Tendencia",
  par.settings = RdBuTheme())
p2 <- levelplot(
  mk[[2]],names.attr = "p-valor",
  par.settings = RdBuTheme())

# Graficando
grid.arrange(p1,p2, ncol = 2)

Lo anteriormente realizado hizo uso de listas de rasters, sin embargo muchas veces la información no se encuentra en este formato. En ocasiones se tiene archivos en formato netcdf.

La página de WorldClim es una buena opción de descarga. La página Terraclimate es más sencillo

Utilizando datos netCDF

library(ncdf4)
dt <- raster::brick("chillón/agg_terraclimate_aet_1958_CurrentYear_GLOBE.nc")

[1] "vobjtovarid4: error #F: I could not find the requsted var (or dimvar) in the file!"
[1] "var (or dimvar) name: crs"
[1] "file name: C:\\Users\\BryanQuispe\\Downloads\\Tendencia\\chillón\\agg_terraclimate_aet_1958_CurrentYear_GLOBE.nc"

plot(dt[[1]])

mk2 <- MKraster(dt, "both")

[1] "Start MKraster: 2023-07-16 02:20:41.832345"
[1] "Loading parameters"
[1] "Done loading parameters"
[1] "Initiating loop operation"
[1] "Populating raster brick"
[1] "Ending MKraster on 2023-07-16 02:20:46.941698"

plot(mk2)

p3 <- levelplot(
  mk2[[1]], names.attr = "Tendencia",
  par.settings = RdBuTheme())
p4 <- levelplot(
  mk2[[2]],names.attr = "p-valor",
  par.settings = RdBuTheme())

# Graficando
grid.arrange(p3,p4, ncol = 2)