#===================================================================================
# Modelamiento hidrológico utilizando como entrada los datos de precipitación de PISCO
# en el modelo hidrológico GR2M
# Profesor: Ing. Carlos A. Fernández Palomino
# Email:cafpxl@gmail.com
#===================================================================================
#===================================================================================
# INSTALAR Y CARGAR PAQUETES NECESARIOS
#===================================================================================
# Instalar paquetes
#install.packages("gstat")
#install.packages("sp")
#install.packages("rgdal")
#install.packages("rgeos")
#install.packages("maptools")
#install.packages("xts")
#install.packages("lubridate")
#install.packages("raster")
#install.packages("corrplot")
#install.packages("cutoffR")
#install.packages("airGR")
#cargar paquetes
library(gstat)
library(sp)
library(rgdal)
## rgdal: version: 1.2-7, (SVN revision 660)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: C:/Program Files/R/R-3.4.0/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: C:/Program Files/R/R-3.4.0/library/rgdal/proj
## Linking to sp version: 1.2-4
library(rgeos)
## rgeos version: 0.3-23, (SVN revision 546)
## GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084
## Linking to sp version: 1.2-4
## Polygon checking: TRUE
library(maptools) # autoloads sp
## Checking rgeos availability: TRUE
library(xts)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
require(lubridate)
## Loading required package: lubridate
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
library(raster)
library(corrplot)
library(cutoffR)
library(airGR)
## Warning: package 'airGR' was built under R version 3.4.2
library(lattice)
library(latticeExtra)
## Loading required package: RColorBrewer
#===================================================================================
# GUARDANDO EL PROYECTO EN EL DIRECTORIO DE TRABAJO Y CREANDO VARIABLES PRINCIPALES
#===================================================================================
#WF <- tk_choose.dir(getwd(), "Seleccionar carpeta")#Carpeta de trabajo lo más cerca de MIS DOCUMENTOS
WF <- "E:/SENAMHI_2017/Curso_R_Hidrologia/Modelamiento_Hidrologico"#Carpeta de trabajo lo más cerca de MIS DOCUMENTOS
setwd(WF)
#===================================================================================
# CARGAR FUNCIONES AUXILIARES
#===================================================================================
source("Script.auxiliares/idw.utilidades.R")
source("Script.auxiliares/utilidades.R")
#===================================================================================
#DEFINIENDO SISTEMA DE COORDENADAS Y VARIABLE A INTERPOLAR
#===================================================================================
utm<-"+init=epsg:32718"# EPSG:32718: WGS 84 / UTM zone 18S #http://spatialreference.org/ref/epsg/32718/
wgs<-"+proj=longlat +datum=WGS84"
#variable = "TEMP"#variable a interpolar
id.comun = "CODIGO" #identificador de estaciones
#===================================================================================
# IMPORTAR Y DAR FORMATO A LOS DATOS;
#===================================================================================
#importando datos diarios
#cambiar 'prec' por 'tmax' o 'tmin'
data.diaria <- as.xts(read.zoo("Datos/prec.csv",header = TRUE,sep = ",",check.names = FALSE))#dataframe de fechas y registros de estaciones
plot.zoo(as.zoo(data.diaria))
#agregando a paso mensual
#data.mensual <- apply.monthly(data.diaria,mean)#por defecto
#FUN=SUM para 'prec' y para 'tmax' y 'tmin' cambiar FUN por MEAN
data.mensual <- diario2mensual_xts_df(data.diaria, FUN=SUM)# tiene en cuenta datos cantidad de datos perdidos
head(data.mensual)
## 679 683 844 809 690 687 812 759
## 2000-01-31 395.4 172.0 159.8 149.4 169.3 NA 119.2 110.4
## 2000-02-29 339.5 92.9 105.4 77.0 110.5 206.4 184.2 180.0
## 2000-03-31 265.1 87.3 58.7 81.9 107.7 NA 112.0 121.9
## 2000-04-30 207.2 NA 2.2 2.9 20.1 39.0 30.6 15.6
## 2000-05-31 75.3 6.4 9.1 1.8 2.0 7.3 7.5 7.2
## 2000-06-30 132.2 8.8 3.7 11.6 1.8 8.7 14.8 6.0
# DEFINIR PERIODO DE ANALISIS AQUI
obs_d <- window(data.mensual[,],start="2000-01-01",end="2015-12-31") #>>>>>>>>>> definir periodo de analisis <<<<<<<<<<<
plot.zoo(as.zoo(obs_d),main="Series de tiempo sin completar")
nonNAs <- function(x) { as.vector(apply(x, 1, function(x) length(which(!is.na(x)))))}
noNAs <- nonNAs(obs_d)#estaciones por mes
#plot para ver los datos observados por mes
plot(cbind(obs_d[,1],noNAs)[,2],type="o",main="Número de datos por dÃa")
#importando datos de estaciones
(data.estacion <- read.csv("Datos/estaciones.csv"))#dataframe de caracteristicas de estaciones
## CODIGO NOMBRE LAT LONG ELEVATION
## 1 679 679 -13.16694 -72.54583 2108
## 2 683 683 -13.31056 -72.12389 3113
## 3 844 844 -13.41611 -71.84972 3222
## 4 809 809 -13.60028 -71.70028 3069
## 5 690 690 -13.61000 -71.56028 3688
## 6 687 687 -13.91694 -71.68361 3227
## 7 812 812 -14.02806 -71.57278 3686
## 8 759 759 -14.25361 -71.23722 3559
str(data.estacion)
## 'data.frame': 8 obs. of 5 variables:
## $ CODIGO : int 679 683 844 809 690 687 812 759
## $ NOMBRE : int 679 683 844 809 690 687 812 759
## $ LAT : num -13.2 -13.3 -13.4 -13.6 -13.6 ...
## $ LONG : num -72.5 -72.1 -71.8 -71.7 -71.6 ...
## $ ELEVATION: int 2108 3113 3222 3069 3688 3227 3686 3559
#importando shapefile de la cuenca
cuenca.shp <- readShapePoly('Datos/Shapefiles/Basin.shp')
## Warning: use rgdal::readOGR or sf::st_read
proj4string(cuenca.shp) <- CRS(utm)
#cuenca.shp<- spTransform(cuenca.shp, CRS(wgs))#transformando projection
plot(cuenca.shp)
save(cuenca.shp,file = "Salidas.R/cuenca.shp.Rdata")# "Salidas.R/tmax.Rdata"
#Calculando el centroide de la cuenca
coordinates(spTransform(cuenca.shp, CRS(wgs)))# anotar los centroies
## [,1] [,2]
## 0 -71.5357 -13.82739
cuenca.centroide <- SpatialPointsDataFrame(coords=coordinates(cuenca.shp), data=data.frame(coordinates(cuenca.shp)),
proj4string=CRS(utm))
plot(cuenca.centroide, col="blue",pch=13,cex=2,add=TRUE)
# crear una grilla vacÃa
raster.base <- raster(extent(cuenca.shp))#baseraster <- raster(extent(-73,-69.9,-15.1,-12.5))
res(raster.base) <- 10000
projection(raster.base) <- CRS(utm)
raster.base[] <- 1
raster.base.sp <- as(raster.base, 'SpatialGridDataFrame')#transformando a objeto sp
# distribucion espacial de las estaciones
dist.esp.est <- data.estacion
coordinates(dist.esp.est) <- ~ LONG + LAT
proj4string(dist.esp.est)=CRS(wgs)
dist.esp.est <- spTransform(dist.esp.est, CRS(utm))
# plot de datos espaciales de la cuenca
plot(raster.base)
plot(cuenca.shp,add=TRUE)
#plot(raster.base.sp,add=TRUE)
plot(dist.esp.est,col="red",add=TRUE)
#===================================================================================
# COMPLETACION,INTERPOLACION Y PROMEDIO AREAL MEDIANTE IDW
# Como output seleccionar "idw.comp", "idw.mapa" o "idw.prom.areal"
#===================================================================================
# Completando datos mediante IDW
idw.comp <- idw.utilidades(obs_d=obs_d, data.estacion=data.estacion,id.comun = "CODIGO", raster.base.sp=raster.base.sp, cuenca.shp=cuenca.shp, utm=utm, wgs=wgs,output = "idw.comp")
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
plot.zoo(as.zoo(idw.comp),main="Completacion con IDW")
# Completando datos mediante cutoffR
# Revisar: Feng, L., Nowak, G., O'Neill, T.J., & Welsh, A.H. (2014). CUTOFF: A spatio-temporal imputation method. Journal of Hydrology, 519(Part D), 3591-3605.
# Correlacion cruzada
library(RColorBrewer)
corrplot(cor(as.matrix(obs_d),use="complete"), method="color", col=colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))(200),
type="upper", order="hclust",
addCoef.col = "black", # Add coefficient of correlation
tl.col="black", tl.srt=45, #Text label color and rotation
# Combine with significance
p.mat = cor.mtest(as.matrix(obs_d)), sig.level = 0.01, insig = "blank",
# hide correlation coefficient on the principal diagonal
diag=FALSE
)
data.cutoff <- data.frame(obs_d,date=index(obs_d),check.names = FALSE)
#data.cutoff.comp <- cutoff(data = data.cutoff)
data.cutoff.comp <- cutoff(data = data.cutoff,method = c("correlation"),corr = "spearman",cutoff = 0.7)# para prec:cutoff = 0.7, tmin:cutoff = 0.8, tmax:cutoff = 0.7
data.mensual.comp <- as.xts(data.cutoff.comp,order.by = index(obs_d))
plot.zoo(as.zoo(data.mensual.comp),main="Completacion con cutoff")
# comparacion de tecnicas de completacion
xyplot(idw.comp, col="red")+# IDW
xyplot(data.mensual.comp, col="blue")+# CUTOFF
xyplot(obs_d, col="black")# OBS
# # Mapas
# idw.mapa <- idw.utilidades(obs_d=data.mensual.comp, data.estacion=data.estacion,id.comun = "CODIGO", raster.base.sp=raster.base.sp, cuenca.shp=cuenca.shp, utm=utm, wgs=wgs,output = "idw.mapa")
# plot(idw.mapa[[1]])
# plot(cuenca.shp,add=TRUE)
# plot(dist.esp.est,col="red",add=TRUE)
# Promedio areal
idw.prom.areal <- idw.utilidades(obs_d=data.mensual.comp, data.estacion=data.estacion,id.comun = "CODIGO", raster.base.sp=raster.base.sp, cuenca.shp=cuenca.shp, utm=utm, wgs=wgs,output = "idw.prom.areal")
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
plot.zoo(as.zoo(idw.prom.areal))
# Guardar Datos
prec <- idw.prom.areal
save(prec,file = "Salidas.R/prec.Rdata")
# Repetir todos los procesos anteriores para la Tmax
# tmin <- idw.prom.areal#tmin or tmax
# save(tmin,file = "Salidas.R/tmin.Rdata")# "Salidas.R/tmax.Rdata"
#
# tmax <- idw.prom.areal#tmin or tmax
# save(tmax,file = "Salidas.R/tmax.Rdata")# "Salidas.R/tmax.Rdata"
#remover todos los objetos y luego repetir para tmin
rm(list = ls())
#===================================================================================
# ESTIMACION DE EVAPOTRANSPIRACION POTENCIAL: cargar la funcion et.hargreaves
#===================================================================================
# Otros metodos revizar el paquete "Evapotranspiration
source("Script.auxiliares/utilidades.R")
load(file = "Salidas.R/tmax.Rdata")
load(file = "Salidas.R/tmin.Rdata")
ET.har <- et.hargreaves(tmin,tmax,lat=-13.82)
plot(ET.har)
save(ET.har,file = "Salidas.R/ETP.har.Rdata")# "Salidas.R/tmax.Rdata"
rm(list = ls())
#===================================================================================
# ESTIMACION DE PRECIPITACION PROMEDIO AREAL PARA LA CUENCA EN BASE AL
# PRODUCTO PISCO - SENAMHI
#===================================================================================
# library(ncdf4)
# utm<-"+init=epsg:32718"# EPSG:32718: WGS 84 / UTM zone 18S #http://spatialreference.org/ref/epsg/32718/
# wgs<-"+proj=longlat +datum=WGS84"
#
# Pisco.prec <- list.files("Datos/PISCO.PREC.MENSUAL",pattern = "*.tif",full.names = TRUE)# lista de geotiffs
# #Pisco.prec.stack <- stack("E:/SENAMHI 2016/BASE_DATOS/PISCO/PISCO_cesar/PiscoMENSUALfinal19812015.nc")# apila los geotiffs
# Pisco.prec.stack <- stack(Pisco.prec)
# plot(Pisco.prec.stack[[1:12]])
# #importando shapefile de la cuenca
# cuenca.shp <- readShapePoly('Datos/Shapefiles/Basin.shp')
# proj4string(cuenca.shp) <- CRS(utm)
# cuenca.shp.wgs <- spTransform(cuenca.shp, CRS(wgs))#transformando projection
#
# plot(Pisco.prec.stack[[1]])
# plot(cuenca.shp.wgs, add=TRUE)
#
# Pisco.prec.prom.cuenca <- extract(Pisco.prec.stack, cuenca.shp.wgs,fun=mean)# extrae el promedio areal para la cuenca
# Pisco.prec.prom.cuenca <- xts(t(Pisco.prec.prom.cuenca), order.by = seq.Date(from = as.Date("1981-01-01"),to=as.Date("2016-12-01"),by="month"))
# plot(Pisco.prec.prom.cuenca)
#
# prec <- window(Pisco.prec.prom.cuenca,start = "2000-01-01",end = "2015-12-31")
# plot(prec,main="Precipitación periodo 2000-2015")
#save(prec,file = "Salidas.R/prec.Rdata")# "Salidas.R/tmax.Rdata"
#load(file = "Salidas.R/prec.Rdata")
#xyplot(prec, col="red")+xyplot(prec2)
#===================================================================================
# CARGAR DATOS DE CAUDALES Y TRANSFORMAR A mm
#===================================================================================
source("Script.auxiliares/utilidades.R")
load(file = "Salidas.R/cuenca.shp.Rdata")
#cargar datos de caudales medios diarios
Q.diario <- as.xts(read.zoo("Datos/Q_diario_vilcanota.csv",header = TRUE,sep = ",",check.names = FALSE))
plot(Q.diario,main="Caudales diarios - estación KM-105")
#Transformar a medios mensuales
Q.mensual <- diario2mensual_xts_df(Q.diario,FUN=MEAN)
plot(Q.mensual, main="Caudales mensuales")
#Tranformar a mm
Q.mensual.mm <- Q.mensual*day(Q.mensual)*24*3600*1000/area(cuenca.shp)
plot(Q.mensual.mm, main="Caudales mensuales [mm]")
Q.mensual.mm <- window(Q.mensual.mm,start = "2000-01-01",end = "2015-12-31")
save(Q.mensual.mm,file = "Salidas.R/Q.mensual.mm.Rdata")# "Salidas.R/tmax.Rdata"
rm(list = ls())
#===================================================================================
# MODELAMIENTO HIDROLOGICO CON GR2M
#===================================================================================
# archivos preprocesados y necesarios
load(file = "Salidas.R/ETP.har.Rdata")
load(file = "Salidas.R/prec.Rdata")
load(file = "Salidas.R/Q.mensual.mm.Rdata")
xyplot(prec, col="red")+# PREC
xyplot(Q.mensual.mm) # Q
## Cargar paquete
library(airGR)
# Datos para GR2M
BasinObs <- data.frame("DatesR"=as.POSIXlt(seq.Date(from = as.Date("2000-01-01"),to=as.Date("2015-12-01"),by="month"),format="%Y-%m-%d"), "P"=as.numeric(prec), "E"=as.numeric(ET.har), "Qmm"=as.numeric(Q.mensual.mm))
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR2M, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
## run period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%m/%Y")=="01/2000"),
which(format(BasinObs$DatesR, format = "%m/%Y")=="12/2010"))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR2M,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## Warning in CreateRunOptions(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel, : Model warm-up period not defined -> default configuration used
## No data were found for model warm-up!
## No warm-up period is used!
## Warning in CreateRunOptions(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel, : Model states initialisation not defined -> default configuration used
## simulation
Param <- c(265.072, 1.040)
OutputsModel <- RunModel_GR2M(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
## results preview
plot_OutputsModel(OutputsModel = OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
## Warning: Deprecated function. Please, use plot.OutputsModel() or plot() on
## an object of class OutputsModel.
## Warning in plot_OutputsModel(OutputsModel = OutputsModel, Qobs = BasinObs
## $Qmm[Ind_Run]): Deprecated "OutputsModel" argument Please, use "x" instead.
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## Crit. NSE[Q] = 0.3400
## preparation of CalibOptions object
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR2M, FUN_CALIB = Calibration_Michel)
## calibration
OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR2M, FUN_CRIT = ErrorCrit_NSE)
## Grid-Screening in progress (
## 0%
## 20%
## 40%
## 60%
## 80%
## 100%)
## Screening completed (9 runs)
## Param = 152.933 , 0.907
## Crit NSE[Q] = 0.3248
## Steepest-descent local search in progress
## Calibration completed (19 iterations, 80 runs)
## Param = 254.678 , 1.485
## Crit NSE[Q] = 0.9248
## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR2M, FUN_CRIT = ErrorCrit_NSE,
FUN_CALIB = Calibration_Michel)
## Grid-Screening in progress (
## 0%
## 20%
## 40%
## 60%
## 80%
## 100%)
## Screening completed (9 runs)
## Param = 152.933 , 0.907
## Crit NSE[Q] = 0.3248
## Steepest-descent local search in progress
## Calibration completed (19 iterations, 80 runs)
## Param = 254.678 , 1.485
## Crit NSE[Q] = 0.9248
## simulation
Param <- OutputsCalib$ParamFinalR #parametros calibrados
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param, FUN = RunModel_GR2M)
## results preview
plot_OutputsModel(OutputsModel = OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
## Warning: Deprecated function. Please, use plot.OutputsModel() or plot() on
## an object of class OutputsModel.
## Warning: Deprecated "OutputsModel" argument Please, use "x" instead.
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## Crit. NSE[Q] = 0.9248
## Validation
## run period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%m/%Y")=="01/2011"),
which(format(BasinObs$DatesR, format = "%m/%Y")=="12/2015"))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR2M,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## Warning in CreateRunOptions(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel, : Model warm-up period not defined -> default configuration used
## The year preceding the run period is used
## Warning in CreateRunOptions(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel, : Model states initialisation not defined -> default configuration used
## simulation
Param <- OutputsCalib$ParamFinalR
OutputsModel <- RunModel_GR2M(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
## results preview
plot_OutputsModel(OutputsModel = OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
## Warning: Deprecated function. Please, use plot.OutputsModel() or plot() on
## an object of class OutputsModel.
## Warning in plot_OutputsModel(OutputsModel = OutputsModel, Qobs = BasinObs
## $Qmm[Ind_Run]): Deprecated "OutputsModel" argument Please, use "x" instead.
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## Crit. NSE[Q] = 0.9615