added 'doReorder_sepYears_sepMembers_combLeadmonths' and related functions

SCRIPTS/doReorder_sepYears_sepMembers_combLeadmonths.R

@@ -6,6 +6,8 @@ library(ecomsUDG.Raccess)
 source(file = "./functions/functionsGeneral.R")
 source(file = "./functions/functionReformat.R")
 source(file = "./functions/functionR2Netcdf.R")
+source(file = "./functions/functionNetcdf2R.R")
+source(file = "./functions/functionConvert.R")
 source(file = "./functions/infoGeneral.R")
 
 submitscript <- FALSE
@@ -16,73 +18,109 @@ if (submitscript) {
   initMonths  <- c(X:X)
   locName    <- 'X'
   resolution<-"0.75"
-  inPath<-sprintf("../DATA/System4_seasonal_15_rev1.0/%s_%s", locName, resolution)
-  outPath<-sprintf("../DATA/System4_seasonal_15_rev1.1REF/%sdeg/%s_noBC", resolution, locName)
+  # inPath<-sprintf("../DATA/System4_seasonal_15_rev1.1REF/%sdeg/%s_noBC", resolution, locName)
+  # outPath<-sprintf("../DATA/System4_seasonal_15_rev1.0/%s_%s", locName, resolution)
 } else {
-  members    <-c(1:2)
-  initYears <-c(1981:1982)
-  initMonths<-1
-  #locName<-"GHA"
-  locName<-"EU"
+  members    <-c(1:15)
+  initYears <-c(1981:2010)
+  initMonths<-c(1:1)
+  locName<-"GHA"
+  #locName<-"EU"
   resolution<-"0.75"
-  inPath<-sprintf("../DATA/System4_seasonal_15_rev1.0/%s_%s", locName, resolution)
-  outPath<-sprintf("../DATA/System4_seasonal_15_rev1.x/%sdeg/%s_noBC", resolution, locName)
+  inPath<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg/%s_noBC", resolution, locName)
+  outPath<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg/%s_noBC_out", resolution, locName)
 }
 
 variables<-names(variableInfo)
 #variables<-c( "pr"    ,  "psl"    , "rsds"   , "rlds" ,   "tasmin","tasmax", "huss" ,   "sfcWind")
+variables<-c( "pr")
 dir.create(outPath, recursive = TRUE, showWarnings = FALSE)
 
 leadMonths<-c(0:6)
 
+sYear<-initYears[1]
+eYear<-initYears[length(initYears)]
+nYears<-length(initYears)
+
+sMember<-members[1]
+eMember<-members[length(members)]
+nMembers<-length(members)
+
+
 print("start")
 for (initMonth in initMonths) {
   for (variableName in variables) {
-    if ('ecomsName' %in% names(variableInfo[[variableName]])) {
-      variableNameECOMS<-variableInfo[[variableName]]$ecomsName
-    } else {
-      variableNameECOMS<-variableName
-    }
-    inFile <- sprintf("%s/%s/%s_forcing_seas15_%s_noBC_E<MEMBERS>_<YEARS>_%02d.RData", inPath, variableNameECOMS, variableNameECOMS, locName, initMonth)
-    print(inFile)
-    RDataAllLeadMonths <- reformat2Bias(members, initYears, initMonth, variableNameECOMS, locName, iFile = inFile)
+    allocatedOutputArray<-FALSE
     
     for (leadMonth in leadMonths) {
+      indexesYearsLeadMonths<-indexesForYearsPerLeadMonth(sYear,eYear,initMonth)
+      
       targetMonth <- getInitTargetInfo( initYear = 0000, initMonth = initMonth, leadMonth = leadMonth )$targetMonth
       targetSYear <- getInitTargetInfo( initYear = initYears[1], initMonth = initMonth, leadMonth = leadMonth )$targetYear
       targetEYear <- getInitTargetInfo( initYear = initYears[length(initYears)], initMonth = initMonth, leadMonth = leadMonth )$targetYear
-      print(sprintf("targetmonth: %s, leadMonth: %s",targetMonth, leadMonth))
+      #   print(sprintf("targetmonth: %s, leadMonth: %s",targetMonth, leadMonth))
       
-      oPrefix <- sprintf("%s/%s_forcing_seas15_%s_noBC_E%02d-%02d_TAR%4d-%4d_%02d_LM%d",
-                         outPath, variableName, locName, 
+      iPrefix <- sprintf("%s/%s_forcing_seas15_%s_noBC_E%02d-%02d_TAR%4d-%4d_%02d_LM%d",
+                         inPath, variableName, locName, 
                          members[1], members[length(members)],
                          targetSYear, targetEYear,
                          targetMonth, leadMonth)
-      print(oPrefix)
+      print(iPrefix)
+      # load(sprintf("%s.RData",iPrefix))
+      RData<-Netcdf2R(sprintf("%s.nc4",iPrefix),variableName)
       
-      RData<-RDataAllLeadMonths[[paste0("LeadMonth_",leadMonth)]]
-      #         print(RData$InitializationDates)
-      #         print(RData$Dates$start)
+      ## Convert Units
+      RData<-convert(RData, toUnit = variableInfo[[variableName]]$units)
       
-      ## set all negative precipitation values 0
-      if (variableName == "pr") {
-        RData$Data[RData$Data < 0] <- 0
+      if (!allocatedOutputArray) {
+        ## allocate arrays
+        finalList<-NULL
+        listNames<-NULL
+        for (iLM in 1:7) {
+          finalList[[iLM]] <- RData
+          listNames<-c(listNames,paste0("LeadMonth_",iLM-1))
+          totalDays<-indexesYearsLeadMonths$eIndexes[as.character(eYear),iLM]
+          
+          finalList[[iLM]]$Data<-array(0,dim=c(nMembers,totalDays,length(RData$xyCoords$y),length(RData$xyCoords$x)))
+          finalList[[iLM]]$Dates$start<-character(length=totalDays)
+          finalList[[iLM]]$Dates$end<-character(length=totalDays)
+        }
+        names(finalList)<-listNames
+        allocatedOutputArray<-TRUE
       }
+      finalList[[paste0("LeadMonth_",leadMonth)]]<-RData
+    }
+    
+    for (initYear in initYears) {
+      indexes<-indexesOfDaysPerMonth(initMonth, 7, initYear = initYear)
       
-      RData$xyCoords$x[]<-round(RData$xyCoords$x[],2)
-      RData$xyCoords$y[]<-round(RData$xyCoords$y[],2)
-      
-      # Check units
-      RData$Variable$varName<-variableName
-      attr(RData$Variable,"standard_name")<-variableInfo[[variableName]]$standardName
-      attr(RData$Variable,"long_name")<-variableInfo[[variableName]]$longName
-      attr(RData$Variable,"units")<-variableInfo[[variableName]]$unitsEcoms
-      
-      ## add some extra attributes
-      attr(RData,"contact") <- "Wietse Franssen (wietse.franssen@wur.nl)"
-      
-      save(RData, file = paste0(oPrefix, ".RData"))
-      #R2Netcdf(paste0(oPrefix, ".nc4"), RData)
+      ## Init structure
+      RData<-finalList[["LeadMonth_0"]]
+      RData$Data<-array(NA,dim=c(indexes[7,"eIndex"],length(RData$xyCoords$y),length(RData$xyCoords$x)))
+      attr(RData$Data,"dimensions")<-c("time", "lat", "lon")
+      RData$Dates$start<-array("",dim=c(indexes[7,"eIndex"]))
+      RData$Dates$end<-array("",dim=c(indexes[7,"eIndex"]))
+
+      for (iMember in members) {
+        ## Fill with data
+        for (iLM in 1:7) {
+          #print(iLM)
+          month<-indexes[iLM,"month"]
+          year<-indexes[iLM,"year"]
+          datesInCurrMonth <- seq(as.Date(as.Date(sprintf("%04d-%02d-01", year, month))),by = "1 day", length=indexes[iLM,"nDaysTotalInMonth"])
+          iDates<-which(is.element(as.Date(finalList[[iLM]]$Dates$start),datesInCurrMonth),TRUE)
+          RData$Dates$start[indexes[iLM,"sIndex"]:indexes[iLM,"eIndex"]]<-finalList[[iLM]]$Dates$start[iDates]
+          RData$Dates$end[indexes[iLM,"sIndex"]:indexes[iLM,"eIndex"]]<-finalList[[iLM]]$Dates$end[iDates]
+          RData$Data[indexes[iLM,"sIndex"]:indexes[iLM,"eIndex"],,]<-finalList[[iLM]]$Data[iMember,iDates,,]
+        }
+        RData$Members<-finalList[[iLM]]$Members[which(initYears==initYear)]
+        RData$InitializationDates<-finalList[[iLM]]$InitializationDates[which(initYears==initYear)]
+
+        ## Save to file...
+        outFile <- sprintf("%s/%s_forcing_seas15_%s_noBC_E%02d_%04d_%02d", outPath, variableName, locName, iMember, initYear, initMonth)
+        #print(outFile)
+        R2Netcdf(outFile = sprintf("%s.nc4",outFile), RData)
+      }
     }
   }
 }

SCRIPTS/functions/functionConvert.R

@@ -9,8 +9,15 @@ convert <- function(RData, fromUnit = NULL, toUnit) {
   }
   if (fromUnit != toUnit) {
     print(paste0("Convert from ", fromUnit, " to ", toUnit))
+    ## START: Exception handling
+    fromUnit_org<-fromUnit
+    toUnit_org<-toUnit
+    if (fromUnit == "kg m-2 s-1") fromUnit <- "mm s-1"
+    if (toUnit_org == "kg m-2 s-1") toUnit <- "mm s-1"
+    ## END: Exception handling
+    
     RData$Data[]<-ud.convert(RData$Data[],fromUnit,toUnit)
-    attr(RData$Variable,"units")<-toUnit
+    attr(RData$Variable,"units")<-toUnit_org
   } else {
     print(paste0("No conversion needed: fromUnit = ", fromUnit, " and toUnit = ", toUnit))
   }

SCRIPTS/functions/functionNetcdf2R.R

@@ -55,7 +55,9 @@ Netcdf2R <- function(inFile, variableName) {
   
   # data
   lData <- ncvar_get(ncid,variableName)
+  
   attr(lData,"dimensions") <- dimensions
+  attr(lData,"dim") <- ncid$var[[variableName]]$size
   
   # xyCoords
   lxyCoords <- list(x = ncid$dim$lon$vals, y = ncid$dim$lat$vals) 
@@ -87,6 +89,8 @@ Netcdf2R <- function(inFile, variableName) {
       attr(RData,names(attributeList[iAttribute])) <- attributeList[[iAttribute]]
     }
   }
+  RData<-FixRDataDimensionorder(RData)
+  
   nc_close(ncid)
   return(RData)
 }
\ No newline at end of file

SCRIPTS/functions/functionsGeneral.R

+## Sets all kinds of RData orders to the fixed order:
+## c("member", "time","lat","lon") or
+## c("time","lat","lon")
+FixRDataDimensionorder<-function(RData) {
+  nDim<-length(dim(RData$Data))
+  
+  order_old<-attr(RData$Data,"dimensions")
+  if (is.null(order_old)) {
+    stop("'dimensions' attribute must be defined!")
+  }
+  
+  if (nDim == 4) {  
+    order_new<-c("member", "time","lat","lon")
+  } else if (nDim == 3) {  
+    order_new<-c("time","lat","lon")
+  } else {
+    stop("number of dimensions should be 3 or 4!")
+  }
+  to<-match(order_new,order_old)
+  dims_old<-dim(RData$Data)
+  RData$Data<-aperm(RData$Data, c(to))
+  attr(RData$Data,"dimensions")<-order_new
+  print(dim(RData$Data))
+  cat("old order: ");print(sprintf("%s (%s)",order_old,dims_old ))
+  cat("new order: ");print(sprintf("%s (%s)",order_new,dim(RData$Data) ))
+  return(RData)
+}
+
 ## Calculation of the index numbers of a set of months
 ## eg: indexes<-indexesOfDaysPerMonth(1,7,1980)
 indexesOfDaysPerMonth <- function(initMonth, nMonths, initYear) {

SCRIPTS/functions/infoGeneral.R

@@ -46,7 +46,7 @@ variableInfo <- list(
   tas     = list(standardName = "air_temperature",     longName = "Near-Surface Air Temperature",               unitsEcoms = "Celsius", units = "K"),
   tasmin  = list(standardName = "air_min_temperature", longName = "Daily Minimum Near-Surface Air Temperature", unitsEcoms = "Celsius", units = "K"),
   tasmax  = list(standardName = "air_max_temperature", longName = "Daily Maximum Near-Surface Air Temperature", unitsEcoms = "Celsius", units = "K"),
-  pr      = list(standardName = "precipitation_flux",  longName = "Precipitation",                              unitsEcoms = "mm day-1", units = "mm s-1", ecomsName = "tp"),
+  pr      = list(standardName = "precipitation_flux",  longName = "Precipitation",                              unitsEcoms = "mm day-1", units = "kg m-2 s-1", ecomsName = "tp"),
   psl     = list(standardName = "air_pressure_at_sea_level", longName = "Sea level pressure",                   unitsEcoms = "Pa", units = "Pa"),
   rsds    = list(standardName = "surface_downwelling_shortwave_flux_in_air", longName = "Surface Downwelling Shortwave Radiation", unitsEcoms = "1/86400^2 W m-2", units = "W m-2"),
   rlds    = list(standardName = "surface_downwelling_longwave_flux_in_air",  longName = "Surface Downwelling Longwave Radiation",  unitsEcoms = "1/86400^2 W m-2", units = "W m-2"),