Última actualización: 03 junio 2020
Se seleccionaron 20 estaciones pluviométricas en el período 2009-2015. Los resultados de la validación cruzada por el inverso de la distancia y por estación se presenta en la siguiente tabla:
xval = readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/lluvia/xval_stat_2009_2015.RDS")
kable(xval,
caption = "Resultados de la validación cruzada por estación") %>%
kable_styling(c("striped", "bordered")) %>%
scroll_box(width = "500px", height = "200px")| LONG | LAT | R2 | bias |
|---|---|---|---|
| 555.1414 | 6192.219 | 0.7483360 | 7.2 |
| 679.2517 | 6240.584 | 0.7974294 | 2.0 |
| 624.5276 | 6225.934 | 0.8983122 | -3.4 |
| 640.1441 | 6221.278 | 0.6936884 | -4.4 |
| 703.2663 | 6196.823 | 0.7309147 | 5.3 |
| 574.7237 | 6227.569 | 0.7997958 | 1.2 |
| 688.7492 | 6207.114 | 0.7420913 | -9.5 |
| 678.7503 | 6166.258 | 0.5853925 | 7.3 |
| 624.6739 | 6237.023 | 0.9090250 | 6.9 |
| 603.3835 | 6232.843 | 0.7511607 | -1.6 |
| 555.2201 | 6205.526 | 0.7914942 | -8.2 |
| 640.7830 | 6202.411 | 0.7601499 | -3.9 |
| 635.9476 | 6250.177 | 0.6974879 | -8.5 |
| 599.4565 | 6210.705 | 0.8100673 | -1.7 |
| 613.3242 | 6214.984 | 0.8630032 | 2.1 |
| 572.7154 | 6206.515 | 0.7911782 | 4.7 |
| 661.8142 | 6194.305 | 0.7848626 | -1.4 |
| 660.9903 | 6199.866 | 0.7713468 | 3.6 |
| 670.1499 | 6248.517 | 0.7781998 | -2.2 |
| 584.5107 | 6187.563 | 0.5573270 | 3.4 |
El rango de valores del coeficiente R2 es de 0.56-0.91 y del sesgo porcentual es de -9.5-7.3.
El código utilizado para la interpolacion por inverso de la distancia, validación cruzada y la puesta en formato SWAT2012 se presenta a continuación.
rm(list = ls())
library(raster)
library(gstat)
# Sistema de coordenadas
utm = CRS("+proj=utm +zone=21 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
utkm = CRS("+proj=utm +zone=21 +south +ellps=WGS84 +datum=WGS84 +units=km +no_defs")
dec = CRS("+init=epsg:4326")
# Coordenadas de estaciones (se dejan 4 estaciones por fuera, producen valores raros en la xval)
xy <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/precip_coord.RDS")
xy = SpatialPoints(xy[3:2], proj4string = dec)
xy = spTransform(xy, utkm)
ccc = !round(data.frame(xy)[,1],4) %in% c(668.5702, 567.6990, 568.8778, 614.0563)
xy = xy[ccc]
p <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/precip_SL_1981_2016.RDS")
p = p[p$date>=as.Date("2009-01-01"),]
p = data.frame(date=p$date, p[,2:29][,ccc])
# Objetivo a interpolar -> Subcuencas nivel 5.
subs1 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/shapefiles_RDS/subs1.RDS")
# Primera parte. Interpolacion y validación cruzada
p_idw = data.frame(date = c(), x=c(), y=c(), var1.pred=c())
idw_cv = data.frame(var1.pred = c(), observed=c(), LONG=c(), LAT=c())
for(t in 1:nrow(p)){
xyz = data.frame(data.frame(xy),
data.frame(z=as.numeric(p[t, 2:ncol(p)])))[,c("LONG","LAT", "z")]
xyz = na.omit(xyz)
xyz = SpatialPointsDataFrame(xyz[,c("LONG","LAT")], data.frame(z = xyz$z),
proj4string = utkm)
if(all(xyz$z==0)){
p_idw = rbind(p_idw, data.frame(date = p$date[t], x=coordinates(subs1)[,1],
y=coordinates(subs1)[,2], var1.pred=0))
}else{
# IDW
p2 = data.frame(krige(z~1,xyz,newdata=subs1, nmax=3))[,c(1,2,3)]
p_idw = rbind(p_idw, data.frame(date=p$date[t], p2))
cv2 = data.frame(krige.cv(z~1,xyz, nmax=3))[,c(1,3,7,8)]
idw_cv = rbind(idw_cv,data.frame(date=p$date[t],cv2))
print(round(t/nrow(p)*100,3))
}
}
saveRDS(idw_cv,"~/R_projects/SWAT_calibration_ETmodis/resultados/idw_cv_2009_2015.RDS")
#Salida gráfica
library(ggplot2)
library(tidyverse)
library(hydroGOF)
xval_stat = idw_cv %>%
group_by(LONG) %>%
mutate(bias=pbias(var1.pred, observed), R2=cor(var1.pred, observed)^2)
xval_stat=unique(xval_stat[,c("LONG", "LAT", "R2", "bias")])
xval_stat %>%
ggplot(aes(x=LONG, y=bias)) + geom_point()
xval_stat %>%
ggplot(aes(x=LONG, y=R2)) + geom_point() +
geom_vline(xintercept = xval_stat$LONG[22])
abline(v=c(xval_stat$LONG[14],xval_stat$LONG[22]))
saveRDS(xval_stat,
"~/R_projects/SWAT_calibration_ETmodis/resultados/lluvia/xval_stat_2009_2015.RDS")
# Archivo pcp
xy_sb = coordinates(subs1)
i=1
p_sub = p_idw[p_idw$x==xy_sb[i,1],c("date","var1.pred")]
for(i in 2:nrow(xy_sb)){
p_sub = cbind(p_sub,p_idw[p_idw$x==xy_sb[i,1],c("var1.pred")])
}
colnames(p_sub) = c("date", paste0("p",1:41))
coordinates(spTransform(subs1,dec))[,2]
Station = paste0("p",1:41)
Lati = coordinates(spTransform(subs1,dec))[,2]
Long = coordinates(spTransform(subs1,dec))[,1]
Elev = rep(50,length(Lati))
library(stringr)
jday = as.numeric(p_sub$date + 1 - as.Date(paste0(format(p_sub$date,"%Y"),"-01-01")))
jday = str_pad(jday,width=3, pad="0")
jday = paste0(format(p_sub$date, "%Y"),jday)
p_sub = p_sub[,-1]
p_sub = round(p_sub,1)
p_sub2 = format(p_sub[,1], digits = 1)
p_sub2 = str_replace_all(p_sub2, " ", "0")
p_sub3 = cbind(p_sub2)
for(i in 2:ncol(p_sub)){
p_sub2 = format(p_sub[,i], digits = 1)
p_sub2 = str_replace_all(p_sub2, " ", "0")
p_sub3 = cbind(p_sub3,p_sub2)
}
p_out = cbind(jday, p_sub3)
Lati = format(round(Lati,1),nsmall=1)
Long = format(round(Long,1),nsmall=1)
p2w = rbind(paste0("Station ", paste(Station, collapse = ","), collapse = ""),
paste0("Lati ", paste(Lati, collapse = ""), collapse = ""),
paste0("Long ", paste(Long, collapse = ""), collapse = ""),
paste0("Elev ", paste(Elev, collapse = " "), collapse = ""),
t(t(apply(p_out,1,function(x){paste0(x, collapse = "")}))))
write.table(p2w,"~/R_projects/SWAT_calibration_ETmodis/resultados/lluvia/p2_2009_2016.txt",
quote = FALSE, sep = "", row.names = FALSE, col.names = FALSE)A work by Rafael Navas