RFmerge CHIRILU in Peru
Carlos Millan
7/19/2021
1. Setting up the environment
2. Loading input data
# Precipitation terrain data
StationsPp <- read.csv.zoo("c_PpStations/baseTotal_01011981_30112020.csv",format="%d/%m/%Y", FUN =as.Date)
StationsGis <- read.csv("c_PpStations/cab_baseTotal981_2018_Final.csv")
basinsShpSf <- read_sf("b_GeoLayers/PacificBasin.shp")
# basinsShp <- shapefile("b_GeoLayers/PacificBasin.shp")
# Check projection (should transform to geography's WGS 84)
basinsShpSf <- sf::st_transform(basinsShpSf,st_crs("+proj=longlat +datum=WGS84 +no_defs"))
# basinsShp <- spTransform(basinsShp, CRS("+init=epsg:4326"))
# Covariates
PISCO10km <- brick("a_PpDB/PISCOO_precip_daily_until_01jul21.nc")## Loading required namespace: ncdf4IMERG10km <- brick("a_PpDB//IMERG_V06B_Early_daily_until_24may20.nc")
PacificDEM <- raster("~/Shapes/HydroshedDEM/HydroSHED.tif") #hydroshed ~90 m
# Domain area
PacificSf.mask <- basinsShpSf[basinsShpSf$NOMBRE %in% c("Cuenca Chillón","Cuenca Rimac","Cuenca LurÃn"),]
# Pacific.mask <- basinsShp[basinsShp$NOMBRE %in% c("Cuenca Chillón","Cuenca Rimac","Cuenca LurÃn"),]
PacificDEM.mask <- crop(PacificDEM,PacificSf.mask)
PacificDEM.mask <- aggregate(PacificDEM.mask,20)3. Basic exploratory data analysis
3.1 Uniformización de la serie
# Filter station by domain area
stations <- StationsGis
stations <- st_as_sf(stations, coords = c('LON', 'LAT'), crs = st_crs("+proj=longlat +datum=WGS84 +no_defs"))
# stations <- st_as_sf(crop(as(stations,"Spatial"),as(PacificSf.mask,"Spatial")))
stations <- st_as_sf(crop(as(stations,"Spatial"),PacificDEM.mask))
plot(as(PacificSf.mask,"Spatial"))
plot(as(stations,"Spatial"),add=T)
StationsPp <- StationsPp[,which(names(StationsPp) %in% stations$COD1)]
# head(stations)
# PISCO.total <- sum(PISCO10km, na.rm= FALSE)/length(unique(lubridate::year(ldates)))
# IMERG.total <- sum(IMERG10km, na.rm= FALSE)/length(unique(lubridate::year(ldates)))
# plot(PISCO.total, main = "PISCO Operativo v2 [2014 - 2019] ", xlab = "Longitude", ylab = "Latitude")
# plot(basinsShp[1],add=T,col="transparent")Preparing input data
Spatial metadata
plot(PacificDEM.mask, main="HydroSHED", xlab="Longitude", ylab="Latitude", col=terrain.colors(255))
plot(PacificSf.mask[1], add=TRUE, col="transparent")
plot(stations[1], add=TRUE, pch = 16, col="black")
6.2 Verification of covariates
nlayers(PISCO10km)## [1] 2191nlayers(IMERG10km)## [1] 2191nrow(StationsPp)## [1] 14579StationsPp <- StationsPp[which(STATIONSdates %in% ldates),]
nrow(StationsPp)## [1] 2191# Homogenize spatial co variables
PISCO10km.mask <- crop(PISCO10km,PacificSf.mask)
IMERG10km.mask <- crop(IMERG10km,PacificSf.mask)
PacificDEM.mask <- crop(PacificDEM.mask,PacificSf.mask)
extent(PISCO10km.mask)## class : Extent
## xmin : -77.2
## xmax : -76
## ymin : -12.3
## ymax : -11.3extent(IMERG10km.mask)## class : Extent
## xmin : -77.2
## xmax : -76
## ymin : -12.3
## ymax : -11.3( extent(PISCO10km.mask) == extent(IMERG10km.mask) )## [1] FALSE( extent(PISCO10km.mask) == extent(PacificDEM.mask) )## [1] FALSE# head(raster::coordinates(PISCO10km.mask))
# head(raster::coordinates(IMERG10km.mask))res(PISCO10km.mask)## [1] 0.1 0.1( res(PISCO10km.mask) == res(IMERG10km.mask) )## [1] TRUE TRUE( res(PISCO10km.mask) == res(PacificDEM.mask) )## [1] FALSE FALSE6.2.1 Optional: Reprojection into another CRS
utm <- CRS("+init=epsg:32718") # WGS 84 / UTM zone 18S
PISCO10km.mask.utm <- projectRaster(from=PISCO10km.mask, crs=utm)
IMERG10km.mask.utm <- projectRaster(from=IMERG10km.mask, crs=utm)
PacificDEM.mask.utm <- projectRaster(from=PacificDEM.mask, crs=utm)stations.utm <- sf::st_transform(stations, crs=32718) # for 'sf' objects# basinsShpSf.utm <- sf::st_transform(basinsShpSf, crs=32718)
# Pacific.mask.utm <- spTransform(Pacific.mask, CRSobj = utm)
PacificSf.mask.utm <- st_transform(PacificSf.mask,crs = 32718)# Get utm station coordinates
st.coords <- st_coordinates(stations.utm)
lon <- st.coords[, "X"]
lat <- st.coords[, "Y"]
StationsPp.utm <- data.frame(ID=stations.utm[["COD1"]], lon=lon, lat=lat)7 Running RFmerge
7.1 Covariates
t1 <- Sys.time()
PISCO10km.mask.utm.res <- resample(PISCO10km.mask.utm,PacificDEM.mask.utm)
t2 <- Sys.time()
IMERG10km.mask.utm.res <- resample(IMERG10km.mask.utm,PacificDEM.mask.utm)
t3 <- Sys.time()
t2-t1## Time difference of 4.613492 secst3-t2## Time difference of 4.760985 secsres(PISCO10km.mask.utm.res)## [1] 1760 1790( res(PISCO10km.mask.utm.res) == res(IMERG10km.mask.utm.res) )## [1] TRUE TRUE( res(PISCO10km.mask.utm.res) == res(PacificDEM.mask.utm) )## [1] TRUE TRUEPISCO10km.mask.utm.res <- crop(PISCO10km.mask.utm.res,PacificSf.mask.utm)
IMERG10km.mask.utm.res <- crop(IMERG10km.mask.utm.res,PacificSf.mask.utm)
PacificDEM.mask.utm <- crop(PacificDEM.mask.utm,PacificSf.mask.utm)
covariates.utm <- list(pisco=PISCO10km.mask.utm.res, imerg=IMERG10km.mask.utm.res,
dem=PacificDEM.mask.utm)
compareRaster(covariates.utm)## [1] TRUEplot(PacificDEM.mask.utm,colNA="red")
plot(PacificSf.mask.utm,add=T, col="transparent")
plot(stations.utm[1],add=T, col="red")
spplot(brick(PISCO10km.mask.utm.res[[1]],IMERG10km.mask.utm.res[[1]]))
#7.2 Setup
# Mask
m <- extent(PacificSf.mask.utm)
m.def = st_sf(ID="extent",st_sfc(st_polygon(list(cbind(c(m[1],m[2],m[2],m[1],m[1]),c(m[3],m[3],m[4],m[4],m[3]))))))
par.nnodes <- min(parallel::detectCores()-1, 45)
drty.out <- file.path("d_Output", "Test3.par")
rfmep <- RFmerge(x=StationsPp, metadata= StationsPp.utm, cov=covariates.utm,
id="ID", lat="lat", lon="lon", mask=m.def, #PacificSf.mask.utm[1]
training=0.8, write2disk=TRUE, drty.out=drty.out,
parallel="parallel", par.nnodes=par.nnodes)## [ Creating the training (80%) and evaluation (20%) datasets ... ]## [ Computing the Euclidean distances to each observation of the training set ...]## ## [ Parallel initialization ... ]## [ Number of cores/nodes detected: 48 ]## [ Number of cores/nodes used : 45 ]## [ Parallelisation finished ! ]ts.path <- paste0(drty.out, "/Ground_based_data/Evaluation/Evaluation_ts.txt")
metadata.path <- paste0(drty.out, "/Ground_based_data/Evaluation/Evaluation_metadata.txt")eval.ts <- read.zoo(ts.path, header = TRUE)
eval.gis <- read.csv(metadata.path)# promoting 'eval.gis' into a spatial object (in order to be plotted)
( eval.gis.utm <- st_as_sf(eval.gis, coords = c('lon', 'lat'), crs = 32719) )## Simple feature collection with 5 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 303536.3 ymin: 8690119 xmax: 362730.5 ymax: 8744682
## Projected CRS: WGS 84 / UTM zone 19S
## ID geometry
## 1 X111083 POINT (303536.3 8744682)
## 2 X111026 POINT (323621.4 8731895)
## 3 X111057 POINT (314626.7 8715249)
## 4 X111058 POINT (333154.1 8690119)
## 5 X111061 POINT (362730.5 8702501)plot(rfmep[[11]], main = "RF-MEP precipitation for 2014-01-11", xlab="Longitude", ylab="Latitude")
plot(PacificSf.mask.utm[1], add=TRUE, col="transparent")
plot(stations.utm[1], add=TRUE, pch = 16, col="blue")
plot(eval.gis.utm, add=TRUE, pch = 16, col="red")
rfmep.total <- sum(rfmep, na.rm= FALSE)
plot(rfmep.total, main = "RF-MEP [2014 - 2019]", xlab = "Longitude", ylab = "Latitude")
plot(PacificSf.mask.utm[1], add=TRUE, col="transparent")
coordinates(eval.gis) <- ~ lon + lat
eval.gis <- eval.gis
rfmep.ts <- t(raster::extract(rfmep, eval.gis))
pisco.ts <- t(raster::extract(PISCO10km.mask.utm, eval.gis))
imerg.ts <- t(raster::extract(IMERG10km.mask.utm, eval.gis))sres <- list(pisco.ts, imerg.ts, rfmep.ts)
nsres <- length(sres)
nstations <- ncol(eval.ts) # -1 pues una estacion se encuentra fuera del area interpolada
tmp <- rep(NA, nstations)
kge.table <- data.frame(ID=eval.gis[["ID"]], PISCO_Op=tmp, IMERG_GPM=tmp, RF_MEP=tmp)
# Computing the KGE between the observed rainfall measured in each one of the raingauges
# of the training dataset and CHIRPSv2, PERSIANN-CDR, the merged product `rfmep`:
for (i in 1:nsres) {
ldates <- time(eval.ts)
lsim <- zoo(sres[[i]], ldates)
kge.table[, (i+1)] = hydroGOF::KGE(sim= lsim, obs= eval.ts, method="2012") # -1 pues una estacion se encuentra fuera del area interpolada
} # FOR end# Boxplot with a graphical comparison
sres.cols <- c("powderblue", "palegoldenrod", "mediumseagreen")
boxplot(kge.table[,2:4], main = "KGE evaluation for Jan - Aug 1983",
xlab = "P products", ylab = "KGE'", ylim = c(-0.5, 1), # horizontal=TRUE,
col=sres.cols)
legend("topleft", legend=c("PISCO_Op", "IMERG_GPM", "RF-MEP"), col=sres.cols,
pch=15, cex=0.8, bty="n")
grid()