Receta
Lizeth C.D.
2023-07-23
Descripción
El presente documento busca reunir todo lo aprendido en un curso de GIS que usa R. Se muestra una primera parte las pruebas de tendencia espacial numérica de los datos de precipitación de estaciones meteorológicas de una zona de estudio comparadas con los datos del producto grillado PISCO. La segunda parte se compone del análisis espacial realizado de la tendencia espacial del producto grillado pisco y 3 indicadores estadísticos en cada estación meteorológica considerada
Análisis de tendencia de series de tiempo
Puntual
Obteniendo data modelada (PISCO), para comparación con estaciones
library(ncdf4)
library(tidyverse)
raster_pp <- raster::brick("PISCOpd.nc")
ce <- read.csv("coords.csv",sep = ";")
names(ce) <- c("estacion","lat","lon","alt")
sp::coordinates(ce) <- ~lon+lat
raster::projection(ce) <- raster::projection(raster_pp)
cells <- raster::extract(
raster_pp[[1]],ce, cellnumbers=T)[,1]
data_cord<-t(raster_pp[cells]) %>%
as.data.frame() %>%
mutate(fecha = seq(
as.Date("1981-01-01"),
as.Date("2016-12-31"),
"days")) %>%
filter(fecha >= as.Date("1981-01-01")&
fecha <= as.Date("2000-12-31"))
names(data_cord) <- c(sort(
ce$estacion,decreasing = T),"fecha")
data_cord <- xts::xts(data_cord[,-7],order.by = data_cord$fecha)Procesando los datos diarios de todas las estaciones escogidas
# creado de un vector de nombres para leer cada excel:
nombres <- list.files(
full.names = FALSE,
pattern = ".xlsx") #%>%
#str_replace("\\.xlsx$","") %>%
#toupper()
# Creado de un vector de nombres para dar nombres a las estaciones
nombres2 <- list.files(
full.names = FALSE,
pattern = ".xlsx") %>%
str_replace("\\.xlsx$","") %>%
toupper()
# Creación de una lista que contenga los datos de todas las estaciones:
data <- lst()
diaria <- lst()
mensual <- lst()
for (i in 1:length(nombres)) {
data[[i]] <- read.xlsx(nombres[i])
names(data[[i]]) <- c("y","d","01",
"02","03","04",
"05","06","07",
"08","09","10",
"11","12")
diaria[[i]] <- pivot_longer(
data = data[[i]],cols = 3:14,
names_to = "m",
values_to = nombres2[i]) %>%
mutate(x=as.Date(paste(
y,m,d,sep = "-"),
format="%Y-%m-%d")) %>%
select(x,nombres2[i]) %>%
arrange(x) %>% na.omit()
}
# Algún día puede servir este código xd
# %>% rename_all(~c("FECHA",nombres))
# Combinando todas las tablas de la lista en una sola:
tabla <- diaria %>% reduce(left_join, by="x")
### Rellenando todos los espacios vacíos:
# creación de una columna de fechas completas
completa <- data.frame(x = seq(min(tabla$x),max(tabla$x),"1 day"))
# agregando las fechas faltantes a la base de datos principal:
tabla <- merge(completa,tabla, by = "x", all.x = TRUE)
# Guardando la tabla:
write.csv(tabla, "estaciones.csv",row.names = F,na = "")Completando datos con KNN (versión python)
import pandas as pd
from sklearn.impute import KNNImputer
# Cargado de la base de datos
df = pd.read_csv("estaciones.csv",
index_col=['x'],
parse_dates=True)
# Estableciendo rango de fechas seleccionado
estaciones = df['01-01-1970':'31-12-2000']
# Estableciendo rango de fechas para el K-Nearest Neighbors
# (vecinos más cercanos)
fechas = pd.date_range('1970-01-01','2000-12-31')
# Método K-Nearest Neighbors (vecinos más cercanos)
imputer = KNNImputer(n_neighbors=6)
df_knn = pd.DataFrame(imputer.fit_transform(
estaciones), columns=df.columns, index=estaciones.index)
estacionesknn = "estaciones_corregido_knn.csv"
df_knn.to_csv("estaciones_corregido_knn.csv", index=True)Preparando las tablas de datos para su comparación
Acondicionando los rangos de fecha:
library(hydroTSM)
library(xts)
library(zoo)
observada <- read.csv("estaciones_corregido_knn.csv")
observada <- observada %>% mutate(x=as.Date(
x,format="%Y-%m-%d")) %>%
filter(x >= as.Date("1981-01-01")&
x <= as.Date("2000-12-31"))
write.csv(observada,"chira_observada.csv", row.names = F)Realizando ajustes finales
library(hydroTSM)
library(xts)
library(zoo)
# Cargando la data observada diaria
observada <- as.xts(
x = read.csv.zoo("chira_observada.csv"))
# Homogenizando los nombres de columnas de lo modelado:
names(data_cord) <- sort(names(observada),
decreasing = T)
modelada <- data_cordConvirtiendo a series mensuales y anuales
# Para la data observada:
observada_m <- hydroTSM::daily2monthly(
observada, FUN=sum, na.rm=TRUE)
observada_y <- hydroTSM::monthly2annual(
observada_m, FUN=sum, na.rm=FALSE)
# Para la data modelada:
modelada_m <- hydroTSM::daily2monthly(
modelada, FUN=sum, na.rm=TRUE)
modelada_y <- hydroTSM::monthly2annual(
modelada_m, FUN=sum, na.rm=FALSE)Realizando test de Mann-kendall
library(Kendall)
library(trend)
# Para la data observada diaria
mk_obsd <- lapply(observada,mk.test)
mk_obsm <- lapply(observada_m,mk.test)
mk_obsy <- lapply(observada_y,mk.test)
# Todas las pruebas resultan en no tendencia:
mk_obsd$AYABACA1##
## Mann-Kendall trend test
##
## data: X[[i]]
## z = 0, n = 7305, p-value = 1
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 0 36155007805 0Realizando la comparación entre la data modelada y observada
library(hydroGOF)
ev_d <- gof(modelada, observada)
ev_m <- gof(modelada_m, observada_m)
ev_y <- gof(modelada_y, observada_y)Generando tablas resumen de los tres indicadores principales PBIAS, RMSE y r.
ind_d <- ev_d[c(4,6,16),]
ind_m <- ev_m[c(4,6,16),]
ind_y <- ev_y[c(4,6,16),]
# Visualizando las tablas creadas:
ind_d## SAUSALCULUCAN1 SAPILLICA1 PANANGA1 MALLARES1 LANCONES1 AYABACA1
## RMSE 7.57 8.88 5.15 6.71 8.15 7.30
## PBIAS % -71.00 2.80 -67.00 -48.60 -36.90 273.90
## r 0.69 0.35 0.63 0.28 0.29 0.64ind_m## SAUSALCULUCAN1 SAPILLICA1 PANANGA1 MALLARES1 LANCONES1 AYABACA1
## RMSE 124.54 110.45 68.89 63.66 117.17 127.61
## PBIAS % -71.00 2.80 -67.00 -48.60 -36.90 273.90
## r 0.93 0.70 0.81 0.89 0.58 0.92ind_y## SAUSALCULUCAN1 SAPILLICA1 PANANGA1 MALLARES1 LANCONES1 AYABACA1
## RMSE 1029.45 683.65 436.10 421.17 729.14 996.73
## PBIAS % -71.00 2.80 -67.00 -48.60 -36.90 273.90
## r 0.88 0.89 0.97 0.94 0.67 0.91# Guardado en 3 hojas de un solo excel:
openxlsx::write.xlsx(lst(ind_d,ind_m,ind_y),
"indicadores.xlsx", rowNames = TRUE)Generación de un shape de las estaciones de monitoreo con la adición de sus indicadores
library(sf)
library(terra)
# Shape para los datos diarios
est_shp <- read.csv("coords.csv",sep = ";")
names(est_shp) <- c("est","lat","lon","alt")
est_shp$est <- names(modelada)
est_shp <- cbind(est_shp,t(ind_d),t(ind_m),t(ind_y))
# Conversión a shp:
est_shp <- est_shp %>% st_as_sf(coords = c("lon","lat"),
crs=4326)Realizando la graficación de la bondad de ajuste para cada estación
ggof(modelada[,1],observada[,6],ftype = "dma",
gofs = c("PBIAS","r","RMSE"),FUN = sum)
ggof(modelada[,2],observada[,5],ftype = "dma",
gofs = c("PBIAS","r","RMSE"),FUN = sum)
ggof(modelada[,3],observada[,4],ftype = "dma",
gofs = c("PBIAS","r","RMSE"),FUN = sum)
ggof(modelada[,4],observada[,3],ftype = "dma",
gofs = c("PBIAS","r","RMSE"),FUN = sum)
ggof(modelada[,5],observada[,2],ftype = "dma",
gofs = c("PBIAS","r","RMSE"),FUN = sum)
ggof(modelada[,6],observada[,1],ftype = "dma",
gofs = c("PBIAS","r","RMSE"),FUN = sum)
Espacial
library(ggspatial)
library(tidyverse)
library(ggrepel)
chira <- st_read("shp/chira.shp")## Reading layer `chira' from data source
## `C:\Users\BryanQuispe\Downloads\Chira_lim\Chira_p\shp\chira.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 21 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -81.15132 ymin: -5.204252 xmax: -79.39046 ymax: -4.192486
## Geodetic CRS: WGS 84ggplot()+
geom_sf(data=chira, fill="turquoise",alpha=0.5)+
geom_sf(data=est_shp,aes(color=`PBIAS %`))+
labs(x="Longitud",y="Latitud",
title="Cuenca Chira",
subtitle = "Estaciones meteorológicas",
caption = "Fuente: Observatorio ANA",
fill="Altitud (m.s.n.m)")+
annotation_scale(location="bl",width_hint=0.2)+
annotation_north_arrow(
location = "tr",height = unit(1,"cm"),
style = north_arrow_fancy_orienteering(
line_width = 0.8))+
theme_minimal()