Deflate

src/Assimilate.cpp

@@ -157,7 +157,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
         {
           Qobs  = _pObservedTS[i]->GetSampledValue(nn);   //mean timestep flow
           Qobs2 = _pObservedTS[i]->GetSampledValue(nn+1); //mean timestep flow
-          
+
           //override initial conditions directly
          if ((nn==1) && (Qobs!=RAV_BLANK_DATA)){
             _pSubBasins[p]->SetQout(Qobs);

src/ChannelXSect.cpp

@@ -389,7 +389,7 @@ double CChannelXSect::GetDiffusivity(const double &Q, const double &SB_slope, co
 {
   ExitGracefullyIf(Q<=0,"CChannelXSect::GetDiffusivity: Invalid channel flowrate",BAD_DATA);
 
-  ///< diffusivity from Roberson et al. 1995, Hydraulic Engineering 
+  ///< diffusivity from Roberson et al. 1995, Hydraulic Engineering
   double slope_mult=1.0;
   double Q_mult    =1.0;
   GetFlowCorrections(SB_slope,SB_n,slope_mult,Q_mult);

src/CustomOutput.cpp

@@ -217,7 +217,57 @@ void CCustomOutput::SetHistogramParams(const double minv,const double maxv, cons
   _hist_max=maxv;
   _nBins=numBins;
 }
-
+//////////////////////////////////////////////////////////////////
+/// \brief Compute the approximate number of time steps for the output array
+/// \param Options [in] Global model options information
+/// \return Number of time steps (int)
+int CCustomOutput::ComputeNumTimeSteps(const optStruct &Options) const
+{
+  // Example assumes Options has julian_start_day, timestep, and custom_interval
+
+  double julian_end_day;
+  int julian_end_year;
+  time_struct start, end;
+
+  JulianConvert(0, Options.julian_start_day, Options.julian_start_year, Options.calendar, start);
+  JulianConvert(Options.duration, Options.julian_start_day, Options.julian_start_year, Options.calendar, end);
+
+  // Compute the end date in Julian day and year format
+  // AddTime(Options.julian_start_day, Options.julian_start_year, Options.duration, Options.calendar, julian_end_day, julian_end_year);
+  // ntime = (int)(ceil((Options.duration+TIME_CORRECTION)/Options.timestep));
+  double dt = Options.timestep;
+  int nsteps = 0;
+  switch(_timeAgg) {
+    case YEARLY:
+      nsteps = int(end.year - start.year);
+      break;
+    case WATER_YEARLY:
+    // wateryr_mo: Starting month of the water year
+      nsteps = int(end.year - start.year);
+      if ((start.month < end.month) & (start.month < Options.wateryr_mo) & (end.month >= Options.wateryr_mo)) nsteps++;
+      if ((start.month > end.month) & (start.month >= Options.wateryr_mo) & (end.month < Options.wateryr_mo)) nsteps--;
+      break;
+    case MONTHLY:
+      nsteps = int((end.year - start.year) * 12 + (end.month - start.month));
+      break;
+    case DAILY:
+      nsteps = int(ceil(Options.duration));
+      break;
+    case EVERY_NDAYS:
+      nsteps = int(ceil(Options.duration / Options.custom_interval));
+      break;
+    case EVERY_TSTEP:
+      nsteps = int(ceil(Options.duration / dt));
+      break;
+    case ENTIRE_SIM:
+      nsteps = 1;
+      break;
+    default:
+      nsteps = 0;
+      break;
+  }
+  return max(0, nsteps);
+}
 ///////////////////////////////////////////////////////////////////
 /// \brief Allocates memory and initialize data storage of a CCustomOutput object
 /// \remarks Called prior to simulation. Determines size of and allocates memory for (member) data[][] array needed in statistical calculations
@@ -554,6 +604,8 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options)
 
   int         retval;                                // error value for NetCDF routines
   size_t      start[1], count[1];                    // determines where and how much will be written to NetCDF
+  size_t      ntime;                                 // Number of time steps
+  size_t      chunksize2[2];
   string      tmp,tmp2,tmp3,tmp4;
 
   bool cant_support=(_aggstat==AGG_RANGE || _aggstat==AGG_95CI || _aggstat==AGG_QUARTILES || _aggstat==AGG_HISTOGRAM);
@@ -574,12 +626,23 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options)
   // time
   // ----------------------------------------------------------
   // (a) Define the DIMENSIONS. NetCDF will hand back an ID for each.
-  retval = nc_def_dim(_netcdf_ID, "time", NC_UNLIMITED, &time_dimid);              HandleNetCDFErrors(retval);
+
+  // TODO: Set dimension size instead of unlimited.
+  ntime = ComputeNumTimeSteps(Options);
+
+  retval = nc_def_dim(_netcdf_ID, "time", ntime, &time_dimid);              HandleNetCDFErrors(retval);
 
   // (b) Define the time variable.
   dimids1[0] = time_dimid;
   retval = nc_def_var(_netcdf_ID, "time", NC_DOUBLE, ndims1,dimids1, &varid_time); HandleNetCDFErrors(retval);
 
+  // Enable deflate compression for time variable (shuffle, zlib, deflate_level)
+  retval = nc_def_var_deflate(_netcdf_ID, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+
+  // Set chunksize to len(time)
+  retval = nc_def_var_chunking(_netcdf_ID, varid_time, NC_CHUNKED, &ntime); HandleNetCDFErrors(retval);
+
+
   // (c) Assign units attributes to the netCDF VARIABLES.
   //     --> converts start day into "hours since YYYY-MM-DD HH:MM:SS"
   char  starttime[200]; // start time string in format 'days since YYY-MM-DD HH:MM:SS'
@@ -617,6 +680,9 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options)
     dimids1[0] = ndata_dimid;
     retval = nc_def_var(_netcdf_ID, group_name.c_str(), NC_STRING, ndims1, dimids1, &varid_grps); HandleNetCDFErrors(retval);
 
+    // Enable deflate compression for group variable
+    retval = nc_def_var_deflate(_netcdf_ID, varid_grps, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+
     //(c) set some attributes to variable "HRUID" or "SBID"
     tmp =long_name;
     tmp2="timeseries_id";
@@ -633,6 +699,14 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options)
     dimids2[1] = ndata_dimid;
     retval = nc_def_var(_netcdf_ID, netCDFtag.c_str(), NC_DOUBLE, ndims2, dimids2, &varid_data);    HandleNetCDFErrors(retval);
 
+    // Enable deflate compression for data variable
+    retval = nc_def_var_deflate(_netcdf_ID, varid_data, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+
+    // Set chunksizes for data variable (time, ndata)
+    chunksize2[0] = ntime; //chunk along time dimension
+    chunksize2[1] = max((size_t)1, min(_nData, (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk
+    retval = nc_def_var_chunking(_netcdf_ID, varid_data, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval);
+
     //(f) set some attributes to variable _netCDFtag
     tmp=_timeAggStr+" "+_statStr+" "+_varName+" "+_spaceAggStr;
     tmp2=_varUnits;
@@ -1267,4 +1341,4 @@ CCustomOutput *CCustomOutput::ParseCustomOutputCommand(char *s[MAXINPUTITEMS], c
   }
 
   return pCustom;
-}
+}

src/CustomOutput.h

@@ -126,7 +126,10 @@ class CCustomOutput
 
   void        SetHistogramParams(const double min,const double max, const int numBins);
 
+  int         ComputeNumTimeSteps(const optStruct &Options) const;
+
   spatial_agg GetSpatialAgg() const;
+
   void        InitializeCustomOutput(const optStruct &Options);
 
   void        WriteFileHeader  (const optStruct &Options);

src/Makefile

@@ -26,8 +26,8 @@ LDFLAGS  :=
 CXXFLAGS += -std=c++11 -fPIC
 
 # OPTION 1) include netcdf - uncomment following two commands (assumes netCDF path = /usr/local):
-#CXXFLAGS += -Dnetcdf
-#LDLIBS   += -L/usr/local -lnetcdf
+CXXFLAGS += -Dnetcdf
+LDLIBS   += -L/usr/local -lnetcdf
 
 # OPTION 1b) include netcdf - for newer MacOS with Apple Silicon (use with option 1 also uncommented)
 #CXXFLAGS += -I/opt/homebrew/include

src/ParsePropertyFile.cpp

@@ -38,10 +38,10 @@ void  AddNewSoilClass       (CSoilClass **&pSoilClasses, soil_struct **&parsed_s
                              string *&soiltags, int &num_parsed_soils, int nConstits,
                              const string name, bool isdefault, CModel *pModel);
 void  AddNewLULTClass       (CLandUseClass **&pLandUseClasses, surface_struct **&parsed_lults,
-                             string *&lulttags, int &num_parsed_lults, 
+                             string *&lulttags, int &num_parsed_lults,
                              const string name, bool isdefault, CModel *pModel);
 void  AddNewVegClass        (CVegetationClass **&pVegClasses, veg_struct **&parsed_vegs,
-                             string *&vegtags, int &num_parsed_vegs, 
+                             string *&vegtags, int &num_parsed_vegs,
                              const string name, bool isdefault, CModel *pModel);
 //////////////////////////////////////////////////////////////////
 /// \brief This method parses the class properties .rvp file
@@ -86,7 +86,7 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
   CSoilProfile     **pProfiles=NULL;
 
   int               num_parsed_aqstacks=0;
-  
+
   bool              invalid_index;
   int               num_read;
   int              *indices=NULL;
@@ -518,7 +518,7 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
       invalid_index=ParsePropArray(p,indices,properties,num_read,soiltags,nParamStrings,num_parsed_soils,aAliases,nAliases);
       ExitGracefullyIf(invalid_index,
                        "ParseClassPropertiesFile: Invalid soiltype code in SoilParameterList command",BAD_DATA);
-      
+
       if (Options.noisy){
         for (int j=0;j<nParamStrings-1;j++){cout<<"  "<<aParamStrings[j+1]<<endl;}
       }
@@ -1520,7 +1520,7 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
     pSoilClasses[c]->AutoCalculateSoilProps        (*parsed_soils[c],*parsed_soils[0],pModel->GetTransportModel()->GetNumConstituents());
   }
   for (int c=1;c<num_parsed_lult;c++) {
-    pLUClasses  [c]->AutoCalculateLandUseProps     (*parsed_surf[c], *parsed_surf[0]);  
+    pLUClasses  [c]->AutoCalculateLandUseProps     (*parsed_surf[c], *parsed_surf[0]);
   }
   for (int c=1;c<num_parsed_terrs;c++){
     pTerrClasses[c-1]->AutoCalculateTerrainProps   ( parsed_terrs[c], parsed_terrs[0]);
@@ -1577,7 +1577,7 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
 
   pModel->CheckForChannelXSectsDuplicates(Options);
 
-  
+
   delete p;
 
   return true;
@@ -1648,7 +1648,7 @@ bool ParsePropArray(CParser          *p,           //parser
         tmpproperties[num_read][j-1]=AutoOrDoubleOrAlias(s[j],pAliases,nAliases);
       }
       DeletePropArray(indices,properties,num_read);
-      
+
       properties=tmpproperties;
       indices=tmpind;
       num_read++;
@@ -1808,7 +1808,7 @@ void AddNewSoilClass(CSoilClass **&pSoilClasses,
 /// \brief dynamically adds a new LULT class to the array of parsed LULT classes.
 //
 void  AddNewLULTClass       (CLandUseClass **&pLandUseClasses, surface_struct **&parsed_lults,
-                             string *&lulttags, int &num_parsed_lults, 
+                             string *&lulttags, int &num_parsed_lults,
                              const string name, bool isdefault, CModel *pModel){
   CLandUseClass *pLC;
   pLC = new CLandUseClass(name, pModel);
@@ -1838,7 +1838,7 @@ void  AddNewLULTClass       (CLandUseClass **&pLandUseClasses, surface_struct **
 /// \brief dynamically adds a new veg class to the array of parsed vegetation classes.
 //
 void  AddNewVegClass        (CVegetationClass **&pVegClasses, veg_struct **&parsed_vegs,
-                             string *&vegtags, int &num_parsed_vegs, 
+                             string *&vegtags, int &num_parsed_vegs,
                              const string name, bool isdefault, CModel *pModel){
   CVegetationClass *pVC;
   pVC = new CVegetationClass(name, pModel);

src/RavenInclude.h

@@ -251,6 +251,8 @@ const double  USE_TEMPLATE_VALUE      =-55555.5;
 const double  NOT_NEEDED              =-66666.6;                                ///< arbitrary value indicating that a non-auto parameter is not needed for the current model configuration
 const double  NOT_NEEDED_AUTO         =-77777.7;                                ///< arbitrary value indicating that a autogeneratable parameter is not needed for the current model configuration
 const double  NETCDF_BLANK_VALUE      =-9999.0;                                 ///< NetCDF flag for blank value
+const int     NETCDF_DEFLATE_LEVEL    = 6;                                      ///< NetCDF deflate level (compression 1-9)
+const int     NETCDF_CHUNKSIZE_MB     = 10;                                     ///< NetCDF chunk size in MB (must be >1 for compression to work)
 const double  RAV_BLANK_DATA          =-1.2345;                                 ///< double corresponding to blank/void data item (also used in input files)
 const double  DIRICHLET_TEMP          =-9999.0;                                 ///< dirichlet concentration flag corresponding to air temperature
 const int     FROM_STATION_VAR        =-55;                                     ///< special flag indicating that NetCDF indices should be looked up from station attribute table

src/Reservoir.cpp

@@ -511,17 +511,17 @@ double  CReservoir::GetLakeConvectionCoeff() const { return _lake_convcoeff; }
 //////////////////////////////////////////////////////////////////
 /// \returns current outflow rate [m3/s]
 //
-double  CReservoir::GetOutflowRate       (const bool adjusted) const { 
+double  CReservoir::GetOutflowRate       (const bool adjusted) const {
   if (adjusted){return max(_Qout+_DAadjust,0.0);}
-  else         {return _Qout;} 
+  else         {return _Qout;}
 }
 
 //////////////////////////////////////////////////////////////////
 /// \returns previous outflow rate [m3/s]
 //
-double CReservoir::GetOldOutflowRate     (const bool adjusted) const { 
+double CReservoir::GetOldOutflowRate     (const bool adjusted) const {
   if (adjusted){return max(_Qout_last+_DAadjust_last,0.0);}
-  else         {return _Qout_last;} 
+  else         {return _Qout_last;}
 }
 
 //////////////////////////////////////////////////////////////////

src/Reservoir.h

@@ -99,7 +99,7 @@ class CReservoir
 
   double        _DAadjust;            //< outflow adjustment - used for reporting overriden flows [m3/s]
   double        _DAadjust_last;       //< outflow adjustment [m3/s] from previous time step
-  
+
   int           _dry_timesteps;      //< number of time steps this reservoir dried out  during simulation
 
   //state variables :

src/StandardOutput.cpp

@@ -2007,6 +2007,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
   int         dimids1[ndims1];                       // array which will contain all dimension ids for a variable
   int         ncid, varid_pre;                       // When we create netCDF variables and dimensions, we get back an ID for each one.
   int         time_dimid, varid_time;                // dimension ID (holds number of time steps) and variable ID (holds time values) for time
+  int         ntime;                                 // Number of time steps
   int         nSim, nbasins_dimid, varid_bsim,varid_bsim2;
   //                                                 // # of sub-basins with simulated outflow, dimension ID, and
   //                                                 // variable to write basin IDs for simulated outflows
@@ -2020,6 +2021,9 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
   string      tmpFilename;
   int         p;                                     // loop over all sub-basins
   string      tmp,tmp2,tmp3;
+  size_t      time_chunksize[1];
+
+
 
   // initialize all potential file IDs with -9 == "not existing and hence not opened"
   _HYDRO_ncid    = -9;   // output file ID for Hydrographs.nc         (-9 --> not opened)
@@ -2054,7 +2058,8 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
   // time
   // ----------------------------------------------------------
   // (a) Define the DIMENSIONS. NetCDF will hand back an ID
-  retval = nc_def_dim(_HYDRO_ncid, "time", NC_UNLIMITED, &time_dimid);  HandleNetCDFErrors(retval);
+  ntime = (int)(ceil((Options.duration+TIME_CORRECTION)/Options.timestep));
+  retval = nc_def_dim(_HYDRO_ncid, "time", ntime, &time_dimid);  HandleNetCDFErrors(retval);
 
   /// Define the time variable. Assign units attributes to the netCDF VARIABLES.
   dimids1[0] = time_dimid;
@@ -2063,6 +2068,13 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
   retval = nc_put_att_text(_HYDRO_ncid, varid_time, "calendar",      strlen("gregorian"), "gregorian"); HandleNetCDFErrors(retval);
   retval = nc_put_att_text(_HYDRO_ncid, varid_time, "standard_name", strlen("time"),      "time");      HandleNetCDFErrors(retval);
 
+  // Enable deflate compression for time variable (shuffle, zlib, deflate_level)
+  retval = nc_def_var_deflate(_HYDRO_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+
+  // Set chunksize to the number of time steps
+  size_t chunksize_time = ntime;
+  retval = nc_def_var_chunking(_HYDRO_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval);
+
   // define precipitation variable
   varid_pre= NetCDFAddMetadata(_HYDRO_ncid, time_dimid,"precip","Precipitation","mm d**-1");
 
@@ -2080,7 +2092,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
     // (b) create dimension "nbasins"
     retval = nc_def_dim(_HYDRO_ncid, "nbasins", nSim, &nbasins_dimid);                             HandleNetCDFErrors(retval);
 
-    // (c) create variable  and set attributes for"basin_name"
+    // (c) create variable  and set attributes for "basin_name"
     dimids1[0] = nbasins_dimid;
     retval = nc_def_var(_HYDRO_ncid, "basin_name", NC_STRING, ndims1, dimids1, &varid_bsim);       HandleNetCDFErrors(retval);
     tmp ="ID of sub-basins with simulated outflows";
@@ -3048,6 +3060,7 @@ int NetCDFAddMetadata2D(const int fileid,const int time_dimid,int nbasins_dimid,
   int    retval;
   int    dimids2[2];
   string tmp;
+  size_t chunksize2[2];
 
   static double fill_val[] = {NETCDF_BLANK_VALUE};
   static double miss_val[] = {NETCDF_BLANK_VALUE};
@@ -3057,6 +3070,15 @@ int NetCDFAddMetadata2D(const int fileid,const int time_dimid,int nbasins_dimid,
 
   // (a) create variable
   retval = nc_def_var(fileid,shortname.c_str(),NC_DOUBLE,2,dimids2,&varid); HandleNetCDFErrors(retval);
+  // Set compression
+  retval = nc_def_var_deflate(fileid, varid, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+  // Set time chunksize to number of time steps
+  retval = nc_inq_dimlen(fileid, time_dimid, &chunksize2[0]); HandleNetCDFErrors(retval);
+   // Set nbasins chunksize to number ensuring that chunks have approximately 10 MB of data (assuming double precision)
+  retval = nc_inq_dimlen(fileid, nbasins_dimid, &chunksize2[1]); HandleNetCDFErrors(retval);
+  chunksize2[1] = max((size_t)1, min(chunksize2[1], (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk
+  // Set chunksize
+  retval = nc_def_var_chunking(fileid, varid, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval);
 
   tmp = "basin_name";