Commit f89fad99 authored by Franssen, Wietse's avatar Franssen, Wietse
Browse files

Committed all old files

parent 31241788
rm(list=ls())
.libPaths(c("/home/WUR/frans004/R/x86_64-unknown-linux-gnu-library/3.1","/cm/shared/apps/R/3.1.2/lib64/R/library"))
library(ecomsUDG.Raccess)
library(WFRTools)
library(ggplot2)
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
library(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
## Calculation of the index numbers of the days of all years
## eg: indexesOfDaysYears(2,1981,1990)
indexesOfDaysYears <- function(month,sYear,eYear) {
nYears<-eYear-sYear+1
indexes <-matrix(0,nYears,5)
colnames(indexes) <- c("sIndex","eIndex","nDaysTotalInMonth","year","month")
dates<-seq(as.Date(paste0(sYear,"-",month,"-1")), by = "1 year", length=nYears)
for (iYear in 1:nYears) {
indexes[iYear,"month"]<-as.numeric(format(dates[iYear], "%m"))
indexes[iYear,"year"]<-as.numeric(format(dates[iYear], "%Y"))
currAndNextMonth<-seq(dates[iYear], by = "1 month", length=2)
indexes[iYear,"nDaysTotalInMonth"]<-difftime(currAndNextMonth[2] ,currAndNextMonth[1] , units = c("days"))
if (indexes[iYear,"year"] == sYear){
indexes[iYear,"eIndex"]<-indexes[iYear,3]
} else {
indexes[iYear,"eIndex"]<-indexes[iYear-1,"eIndex"]+indexes[iYear,"nDaysTotalInMonth"]
}
indexes[iYear,"sIndex"]<-indexes[iYear,"eIndex"]-indexes[iYear,"nDaysTotalInMonth"]+1
}
return(indexes)
}
domainName<-c( "GHA", "EU")
lonmin<- c( 27.75, -24.75)
lonmax<- c( 49.25, 39.75)
latmin<- c(-12.25, 33.25)
latmax<- c( 18.25, 71.75)
barLimits = c(0,10) # EU TP
barLimitsDiff = c(-2.5,2.5)
barLimitsDiff = c(-0.2,0.2)
sYear <- 1981
eYear <- 2010
#eYear <- 1982
## 1 2 3 4 5 6 7 8 9
variables <- c("tas", "tasmax","tasmin", "tp", "psl", "rsds", "rlds", "huss", "wss") # psl is not available as WFDEI
variables <- c("tas", "tasmax","tasmin", "pr", "psl", "rsds", "rlds", "huss", "wss") # psl is not available as WFDEI
variableNames <- c("temperature", "maximum temperature","minimum temperature", "precipitation", "psl", "rsds", "rlds", "huss", "wss") # psl is not available as WFDEI
variableUnits <- c("[Celsius]", "[Celsius]", "[Celsius]", "[mm day-1]", "psl", "rsds", "rlds", "huss", "wss") # psl is not available as WFDEI
sMember=1
eMember=15
#eMember=1
#sourcePath <- "../DATA_RAW_R-BIASFORMAT/%s_0.50"
#destPath <- "../DATA_CORR_R-BIASFORMAT/%s_0.50"
sourcePath <- "../DATA_RAW_R/%s_0.50"
destPath <- "../DATA_CORR_R/%s_0.50"
sourcePath <- "/home/WUR/frans004/L_DATA/PROJECT_DATA/EUPORIAS/System4_seasonal_15/0.50deg/%s_noBC"
#destPath <- "/home/WUR/frans004/L_DATA/PROJECT_DATA/EUPORIAS/System4_seasonal_15/0.50deg/%s_BC"
#destPath <- "/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA_CORR_NETCDF/%s_0.50/final"
#destPath <- "/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA_CORR_NETCDF/%s_0.50_01mmDrizzle/final_pr"
#destPath <- "/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA_CORR_NETCDF/%s_0.50_0mm/final_pr"
destPath <- "/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA_CORR_NETCDF/%s_0.50_01mm/final_complete"
#destPath <- "/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA_CORR_NETCDF/%s_0.50/final_pr"
for (iDomain in 2:2) {
#for (var in c(1,2,3,4,6,7,8,9)) {
for (var in c(4)) {
for (iMonth in 12:12) {
#for (iMonth in 1:12) {
for (iLM in 0:6) {
## Calculate StartMonth of 7 months Run
iStartMonth<-iMonth-iLM
iSeason<-iLM+1
if (iStartMonth<1) {iStartMonth = iStartMonth + 12}
print(paste0("******** Curr-Month: ", iMonth, ", LM: ", iLM, ", S4-StartMonth: ", iStartMonth, ", Season (old): ", iSeason))
dayIndexes<-indexesOfDaysYears(iMonth,sYear,eYear)
tLength<-dayIndexes[eYear-sYear+1,2]
print(dayIndexes)
print(paste0("total Tsteps of all seas15 files: ", tLength))
firstFile<-TRUE
for (iMember in sMember:eMember) {
i<-1
for (iYear in sYear:eYear) {
# ii<-(iYear-sYear)+1
start<-dayIndexes[i,1]
end<-dayIndexes[i,2]
print(paste0("startTstep:", start, " ,endTstep:", end))
## Open SYS4 1
nameFileSys4 <- sprintf(paste0(sourcePath,"/forcing_seas15_%s_noBC_E%02d_%4d_%02d.nc4"),
domainName[iDomain],
domainName[iDomain],
iMember,
iYear,
iStartMonth)
print(paste0("Opening SYS4: ", nameFileSys4))
sys4<-ncLoad(nameFileSys4,
timesteps = "all",
var =variables[var])
sys4_noBC<-sys4
monthlistTmp<-as.numeric(strftime(as.Date(sys4_noBC$Dates$start), format="%m"))
# print(monthlistTmp)
sys4_noBC$Dates$start<-sys4_noBC$Dates$start[monthlistTmp==iMonth]
sys4_noBC$Dates$end<-sys4_noBC$Dates$end[monthlistTmp==iMonth]
sys4_noBC$Data<-sys4_noBC$Data[monthlistTmp==iMonth,,]
if (firstFile) {
sys4_noBCfinal<-sys4_noBC
sys4_noBCfinal$Data<-array(dim = c((eMember-sMember)+1,tLength,dim(sys4_noBC$Data)[2:3]))
}
sys4_noBCfinal$Data[iMember,start:end,,] <-sys4_noBC$Data
rm(sys4)
## Open SYS4 2
nameFileSys4 <- sprintf(paste0(destPath,"/forcing_seas15_%s_BC_E%02d_%4d_%02d.nc4"),
domainName[iDomain],
domainName[iDomain],
iMember,
iYear,
iStartMonth)
# nameFileSys4 <- sprintf(paste0(destPath,"/%s_forcing_seas15_%s_ref__E%02d_%4d_%02d.nc4"),
# domainName[iDomain],
# variables[var],
# domainName[iDomain],
# iMember,
# iYear,
# iStartMonth)
print(paste0("Opening SYS4: ", nameFileSys4))
sys4<-ncLoad(nameFileSys4,
timesteps = "all",
var =variables[var])
sys4_BC<-sys4
monthlistTmp<-as.numeric(strftime(as.Date(sys4_BC$Dates$start), format="%m"))
sys4_BC$Dates$start<-sys4_BC$Dates$start[monthlistTmp==iMonth]
sys4_BC$Dates$end<-sys4_BC$Dates$end[monthlistTmp==iMonth]
sys4_BC$Data<-sys4_BC$Data[monthlistTmp==iMonth,,]
if (firstFile) {
sys4_BCfinal<-sys4_BC
sys4_BCfinal$Data<-array(dim = c((eMember-sMember)+1,tLength,dim(sys4_BC$Data)[2:3]))
}
sys4_BCfinal$Data[iMember,start:end,,] <-sys4_BC$Data
firstFile<-FALSE
i<-i+1
}
}
sys4_noBC<-sys4_noBCfinal
sys4_BC<-sys4_BCfinal
## Open OBS
nameFileObs <- sprintf("../DATA_RAW_R/%s_WFD/%s_wfdei_%s_%s_%02d.Rdata",
domainName[iDomain],
variables[var],
1981,
2010,
iMonth)
print(paste0("Opening OBS: ", nameFileObs))
load(nameFileObs);
## regrid obs
print("Regridding...")
obs_new<-obs
newGrid<-getGrid(sys4)
indexesLat<-c(which(obs$xyCoords$y==newGrid$y[1]):which(obs$xyCoords$y==newGrid$y[2]))
indexesLon<-c(which(obs$xyCoords$x==newGrid$x[1]):which(obs$xyCoords$x==newGrid$x[2]))
obs_new$xyCoords<-sys4$xyCoords
obs_new$Data<-obs$Data[,indexesLat,indexesLon]
attr(obs_new$Data,"dimensions") <- c("time","lat","lon")
obs<-obs_new
#}
## Conversion:
if (variables[var] == "tas" || variables[var] == "tasmin" || variables[var] == "tasmax" ) {
sys4_BC$Data <- sys4_BC$Data-273.15
sys4_noBC$Data <- sys4_noBC$Data-273.15
barLimits = c(-10,10)
}
## Conversion:
if (variables[var] == "pr") {
sys4_BC$Data <- sys4_BC$Data*86400
sys4_noBC$Data <- sys4_noBC$Data*86400
barLimits = c(0,10) # EU TP
}
## Check-Correction:
print("obs_mean...")
obs_mean<-obs;obs_mean$Data<-apply(obs$Data, c(2,3), mean)
# p_obs<-ggplotje(obs_mean, title = paste0("[", variables[var], "] Observed\nmean of month: ", iMonth))
p_obs<-ggplotje(obs_mean, title = paste0(variableNames[var], " ", variableUnits[var], "\nWFDEI (", month.name[iMonth], ")"),barLimits = barLimits)
print("BC_mean...")
sys4_BC_mean<-sys4_BC;sys4_BC_mean$Data<-apply(sys4_BC$Data, c(1,3,4), mean)
p_BC<-ggplotje(sys4_BC_mean, title = paste0(variableNames[var], " ", variableUnits[var], "\nseas15_BC (", month.name[iMonth], ", LM:", iLM, "/SM:", iStartMonth, "/S:", iSeason, ")"),barLimits = barLimits)
print("noBC_mean...")
sys4_noBC_mean<-sys4_noBC;sys4_noBC_mean$Data<-apply(sys4_noBC$Data, c(1,3,4), mean)
p_noBC<-ggplotje(sys4_noBC_mean, title = paste0(variableNames[var], " ", variableUnits[var], "\nseas15_noBC (", month.name[iMonth], ", LM:", iLM, "/SM:", iStartMonth, "/S:", iSeason, ")"),barLimits = barLimits)
sys4_BC_mean2<-sys4_BC
sys4_BC_mean2$Data<-apply(sys4_BC_mean$Data, c(2,3), mean)
obsdiff<-obs_mean
obsdiff$Data<-obs_mean$Data-sys4_BC_mean2$Data
# p_obsVsBC<-ggplotje(obsdiff, title = "Observed vs BC")
p_obsVsBC<-ggplotje(obsdiff, title = paste0(variableNames[var], " ", variableUnits[var], "\nWFDEI - seas15_BC"),barLimits = barLimitsDiff,barDiff=TRUE,backgroudColor = "grey85")
sys4_noBC_mean2<-sys4_BC
sys4_noBC_mean2$Data<-apply(sys4_noBC_mean$Data, c(2,3), mean)
obsdiffnoBC<-obs_mean
obsdiffnoBC$Data<-obs_mean$Data-sys4_noBC_mean2$Data
# p_obsVsnoBC<-ggplotje(obsdiffnoBC, title = "Observed vs noBC")
p_obsVsnoBC<-ggplotje(obsdiffnoBC, title = paste0(variableNames[var], " ", variableUnits[var], "\nWFDEI - seas15_noBC"),barLimits = barLimitsDiff,barDiff=TRUE,backgroudColor = "grey85")
densPlot <-function(lon = NULL, lat = NULL) {
data_obs<-obs$Data[,lat==obs$xyCoords$y,lon==obs$xyCoords$x]
data_BC<-sys4_BC$Data[sMember:eMember,,lat==obs$xyCoords$y,lon==obs$xyCoords$x]
data_noBC<-sys4_noBC$Data[sMember:eMember,,lat==obs$xyCoords$y,lon==obs$xyCoords$x]
dat <- data.frame(legend = factor(c(rep("1: OBS",length(as.vector(data_obs))),
rep("3: noBC",length(as.vector(data_noBC))),
rep("2: BC",length(as.vector(data_BC))))),
data = c(as.vector(data_obs),
as.vector(data_noBC),
as.vector(data_BC)))
nBins<-20
bins<-seq(min(dat$data),max(dat$data),(max(dat$data)-min(dat$data))/(nBins+1))
# Density plots
# pplot<-ggplot(dat, aes(x=data, colour=legend, fill=legend)) +
pplot<-ggplot(dat, aes(x=data, fill=legend)) +
#geom_histogram(binwidth=1, position="dodge",aes(y = ..density..)) +
geom_histogram(breaks=bins, position="dodge",aes(y = ..density..)) +
#geom_density() +
ggtitle(paste0("lon: ",lon, " lat: ", lat))
return(pplot)
}
if (domainName[iDomain] == "GHA") {
p4<-densPlot(30.25,-2.75)
p5<-densPlot(37.25,-2.75)
p6<-densPlot(37.25,10.75)
} else {
# p4<-densPlot(-5.25,40.75)
# p5<-densPlot(30.25,50.75)
# p6<-densPlot(10.25,46.25)
p4<-densPlot(19.25,54.25)
p5<-densPlot(27.25,43.75)
p6<-densPlot(10.25,46.25)
}
nameFilePicture <- sprintf(paste0("./seperateNetCDFFiles_LM/%s_seperateNetCDFFiles_forcing_seas15_%s_E%02d-%02d_%4d-%4d_%02d_LM%1d_01mm_completeXX"),
variables[var],
domainName[iDomain],
sMember,
eMember,
sYear,
eYear,
iMonth,
iLM)
print(paste0("Saving picture : ", nameFilePicture))
png(paste0(nameFilePicture,'.png'), width = 1920, height = 1080)
multiplot(p_obs, p_noBC, p_BC, p_obs, p_obsVsnoBC, p_obsVsBC, p4 ,p5, p6, cols=3)
dev.off()
## SAVE DATA
nameFileDataObs <- sprintf(paste0("./seperateNetCDFFiles_LM_DATA/obs_%s_seperateNetCDFFiles_forcing_seas15_%s_E%02d-%02d_%4d-%4d_%02d_LM%1d_01mm_completeXX.Rdata"),
variables[var],
domainName[iDomain],
sMember,
eMember,
sYear,
eYear,
iMonth,
iLM)
nameFileDataBC <- sprintf(paste0("./seperateNetCDFFiles_LM_DATA/BC_%s_seperateNetCDFFiles_forcing_seas15_%s_E%02d-%02d_%4d-%4d_%02d_LM%1d_01mm_completeXX.Rdata"),
variables[var],
domainName[iDomain],
sMember,
eMember,
sYear,
eYear,
iMonth,
iLM)
nameFileDataNoBC <- sprintf(paste0("./seperateNetCDFFiles_LM_DATA/noBC_%s_seperateNetCDFFiles_forcing_seas15_%s_E%02d-%02d_%4d-%4d_%02d_LM%1d_01mm_completeXX.Rdata"),
variables[var],
domainName[iDomain],
sMember,
eMember,
sYear,
eYear,
iMonth,
iLM)
print(paste0("Saving data : ", nameFileDataObc))
save(data_obs,file=nameFileDataObs)
print(paste0("Saving data : ", nameFileDataBC))
save(data_BC,file=nameFileDataBC)
print(paste0("Saving data : ", nameFileDataNoBC))
save(data_noBC,file=nameFileDataNoBC)
} ## end iMonth
} ## end iSeason
} ## end var
} ## end iDomain
rm(list=ls())
.libPaths(c("/home/WUR/frans004/R/x86_64-unknown-linux-gnu-library/3.1","/cm/shared/apps/R/3.1.2/lib64/R/library"))
library(ecomsUDG.Raccess)
library(WFRTools)
library(ggplot2)
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
library(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
## Calculation of the index numbers of the days of all years
## eg: indexesOfDaysYears(2,1981,1990)
indexesOfDaysYears <- function(month,sYear,eYear) {
nYears<-eYear-sYear+1
indexes <-matrix(0,nYears,5)
colnames(indexes) <- c("sIndex","eIndex","nDaysTotalInMonth","year","month")
dates<-seq(as.Date(paste0(sYear,"-",month,"-1")), by = "1 year", length=nYears)
for (iYear in 1:nYears) {
indexes[iYear,"month"]<-as.numeric(format(dates[iYear], "%m"))
indexes[iYear,"year"]<-as.numeric(format(dates[iYear], "%Y"))
currAndNextMonth<-seq(dates[iYear], by = "1 month", length=2)
indexes[iYear,"nDaysTotalInMonth"]<-difftime(currAndNextMonth[2] ,currAndNextMonth[1] , units = c("days"))
if (indexes[iYear,"year"] == sYear){
indexes[iYear,"eIndex"]<-indexes[iYear,3]
} else {
indexes[iYear,"eIndex"]<-indexes[iYear-1,"eIndex"]+indexes[iYear,"nDaysTotalInMonth"]
}
indexes[iYear,"sIndex"]<-indexes[iYear,"eIndex"]-indexes[iYear,"nDaysTotalInMonth"]+1
}
return(indexes)
}
domainName<-c( "GHA", "EU")
lonmin<- c( 27.75, -24.75)
lonmax<- c( 49.25, 39.75)
latmin<- c(-12.25, 33.25)
latmax<- c( 18.25, 71.75)
barLimits = c(0,10) # EU TP
barLimitsDiff = c(-2.5,2.5)
barLimitsDiff = c(-0.2,0.2)
sYear <- 1981
eYear <- 2010
#eYear <- 1982
## 1 2 3 4 5 6 7 8 9
variables <- c("tas", "tasmax","tasmin", "tp", "psl", "rsds", "rlds", "huss", "wss") # psl is not available as WFDEI
variables <- c("tas", "tasmax","tasmin", "pr", "psl", "rsds", "rlds", "huss", "wss") # psl is not available as WFDEI
variableNames <- c("temperature", "maximum temperature","minimum temperature", "precipitation", "psl", "rsds", "rlds", "huss", "wss") # psl is not available as WFDEI
variableUnits <- c("[Celsius]", "[Celsius]", "[Celsius]", "[mm day-1]", "psl", "rsds", "rlds", "huss", "wss") # psl is not available as WFDEI
sMember=1
eMember=15
#eMember=1
#sourcePath <- "../DATA_RAW_R-BIASFORMAT/%s_0.50"
#destPath <- "../DATA_CORR_R-BIASFORMAT/%s_0.50"
sourcePath <- "../DATA_RAW_R/%s_0.50"
destPath <- "../DATA_CORR_R/%s_0.50"
sourcePath <- "/home/WUR/frans004/L_DATA/PROJECT_DATA/EUPORIAS/System4_seasonal_15/0.50deg/%s_noBC"
#destPath <- "/home/WUR/frans004/L_DATA/PROJECT_DATA/EUPORIAS/System4_seasonal_15/0.50deg/%s_BC"
#destPath <- "/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA_CORR_NETCDF/%s_0.50/final"
#destPath <- "/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA_CORR_NETCDF/%s_0.50_01mmDrizzle/final_pr"
#destPath <- "/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA_CORR_NETCDF/%s_0.50_0mm/final_pr"
destPath <- "/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA_CORR_NETCDF/%s_0.50_01mm/final_complete"
#destPath <- "/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA_CORR_NETCDF/%s_0.50/final_pr"
for (iDomain in 2:2) {
#for (var in c(1,2,3,4,6,7,8,9)) {
for (var in c(4)) {
for (iMonth in X:X) {
#for (iMonth in 1:12) {
for (iLM in 0:6) {
## Calculate StartMonth of 7 months Run
iStartMonth<-iMonth-iLM
iSeason<-iLM+1
if (iStartMonth<1) {iStartMonth = iStartMonth + 12}
print(paste0("******** Curr-Month: ", iMonth, ", LM: ", iLM, ", S4-StartMonth: ", iStartMonth, ", Season (old): ", iSeason))
dayIndexes<-indexesOfDaysYears(iMonth,sYear,eYear)
tLength<-dayIndexes[eYear-sYear+1,2]
print(dayIndexes)
print(paste0("total Tsteps of all seas15 files: ", tLength))
firstFile<-TRUE
for (iMember in sMember:eMember) {
i<-1
for (iYear in sYear:eYear) {
# ii<-(iYear-sYear)+1
start<-dayIndexes[i,1]
end<-dayIndexes[i,2]
print(paste0("startTstep:", start, " ,endTstep:", end))
## Open SYS4 1
nameFileSys4 <- sprintf(paste0(sourcePath,"/forcing_seas15_%s_noBC_E%02d_%4d_%02d.nc4"),
domainName[iDomain],
domainName[iDomain],
iMember,
iYear,
iStartMonth)
print(paste0("Opening SYS4: ", nameFileSys4))
sys4<-ncLoad(nameFileSys4,
timesteps = "all",
var =variables[var])
sys4_noBC<-sys4
monthlistTmp<-as.numeric(strftime(as.Date(sys4_noBC$Dates$start), format="%m"))
# print(monthlistTmp)
sys4_noBC$Dates$start<-sys4_noBC$Dates$start[monthlistTmp==iMonth]
sys4_noBC$Dates$end<-sys4_noBC$Dates$end[monthlistTmp==iMonth]
sys4_noBC$Data<-sys4_noBC$Data[monthlistTmp==iMonth,,]
if (firstFile) {
sys4_noBCfinal<-sys4_noBC
sys4_noBCfinal$Data<-array(dim = c((eMember-sMember)+1,tLength,dim(sys4_noBC$Data)[2:3]))
}
sys4_noBCfinal$Data[iMember,start:end,,] <-sys4_noBC$Data
rm(sys4)
## Open SYS4 2
nameFileSys4 <- sprintf(paste0(destPath,"/forcing_seas15_%s_BC_E%02d_%4d_%02d.nc4"),
domainName[iDomain],
domainName[iDomain],
iMember,
iYear,
iStartMonth)
# nameFileSys4 <- sprintf(paste0(destPath,"/%s_forcing_seas15_%s_ref__E%02d_%4d_%02d.nc4"),
# domainName[iDomain],
# variables[var],
# domainName[iDomain],
# iMember,
# iYear,
# iStartMonth)
print(paste0("Opening SYS4: ", nameFileSys4))
sys4<-ncLoad(nameFileSys4,
timesteps = "all",
var =variables[var])
sys4_BC<-sys4
monthlistTmp<-as.numeric(strftime(as.Date(sys4_BC$Dates$start), format="%m"))
sys4_BC$Dates$start<-sys4_BC$Dates$start[monthlistTmp==iMonth]
sys4_BC$Dates$end<-sys4_BC$Dates$end[monthlistTmp==iMonth]
sys4_BC$Data<-sys4_BC$Data[monthlistTmp==iMonth,,]
if (firstFile) {
sys4_BCfinal<-sys4_BC
sys4_BCfinal$Data<-array(dim = c((eMember-sMember)+1,tLength,dim(sys4_BC$Data)[2:3]))
}
sys4_BCfinal$Data[iMember,start:end,,] <-sys4_BC$Data
firstFile<-FALSE
i<-i+1
}
}
sys4_noBC<-sys4_noBCfinal
sys4_BC<-sys4_BCfinal
## Open OBS
nameFileObs <- sprintf("../DATA_RAW_R/%s_WFD/%s_wfdei_%s_%s_%02d.Rdata",
domainName[iDomain],
variables[var],
1981,
2010,
iMonth)
print(paste0("Opening OBS: ", nameFileObs))
load(nameFileObs);
## regrid obs
print("Regridding...")
obs_new<-obs
newGrid<-getGrid(sys4)
indexesLat<-c(which(obs$xyCoords$y==newGrid$y[1]):which(obs$xyCoords$y==newGrid$y[2]))
indexesLon<-c(which(obs$xyCoords$x==newGrid$x[1]):which(obs$xyCoords$x==newGrid$x[2]))
obs_new$xyCoords<-sys4$xyCoords
obs_new$Data<-obs$Data[,indexesLat,indexesLon]
attr(obs_new$Data,"dimensions") <- c("time","lat","lon")
obs<-obs_new
#}
## Conversion: