Merge pull request #468 from tobeycarman/new-fire-inputs

Flesh out functions in create_region_inputs.py for building explicit fire inputs.

include/Cohort.h

@@ -111,9 +111,12 @@ public :
   void set_state_from_restartdata();
   void set_restartdata_from_state();
 
-  void load_proj_climate(const std::string&); // Overwrites/fills containers with data from projected climate file
-  void load_proj_co2(const std::string& proj_co2_file); // Overwrites/fills co2 container with data from projected file
-
+  // Overwrites/fills climate containers with data from projected climate file
+  void load_proj_climate(const std::string&);
+  // Overwrites/fills co2 container with data from projected co2 input file
+  void load_proj_co2(const std::string& proj_co2_file);
+  // Overwrites/fills explicit fire data containers with data from the projected fire input file.
+  void load_proj_explicit_fire(const std::string& proj_exp_fire_file);
 private:
 
   Integrator vegintegrator[NUM_PFT];

include/WildFire.h

@@ -33,6 +33,8 @@ class WildFire {
 
   ~WildFire();
 
+  void load_projected_explicit_data(const std::string& exp_fname, int y, int x);
+
   int getFRI();
 
   void setCohortData(CohortData* cdp);
@@ -95,7 +97,7 @@ class WildFire {
   std::vector<int> exp_burn_mask;
   std::vector<int> exp_jday_of_burn;
   std::vector<int> exp_fire_severity;
-  std::vector<int> exp_area_of_burn;
+  std::vector<int64_t> exp_area_of_burn;
 
 
 

scripts/input_util.py

@@ -314,6 +314,35 @@ def climate_ts_plot(args):
 
   gs = gridspec.GridSpec(ROWS, COLS)
 
+  def try_infer_time_axis(nc_dataset):
+    '''
+    Tries to build a DatetimeIndex using the 'time' variable 
+    and units of the nc_dataset.
+
+    Parameters
+    ==========
+    nc_dataset : an open netCDF4 dataset. 
+
+    Returns
+    =======
+    idx : a pandas.DatetimeIndex or simple list index starting at 0
+    '''
+
+    try:
+      tV = nc_dataset.variables['time']
+      idx = pd.DatetimeIndex(pd.date_range(
+        start=nc.num2date(tV[0], tV.units, tV.calendar).strftime(),
+        end=nc.num2date(tV[-1], tV.units, tV.calendar).strftime(),
+        freq='MS') # <-- month starts)
+      )
+
+    except RuntimeError as e:
+      print(e)
+      print("Unable to figure out time axis values. Defaulting to basic range index.")
+      idx = np.arange(0, len(nc_dataset))
+
+    return idx
+
   if args.type == 'spatial-temporal-summary':
 
     if args.stitch:
@@ -446,9 +475,14 @@ def climate_ts_plot(args):
       hds = nc.Dataset(os.path.join(args.input_folder, CLIMATE_FILES[0]))
       pds = nc.Dataset(os.path.join(args.input_folder, CLIMATE_FILES[1]))
 
+      h_timeidx = try_infer_time_axis(hds)
+      p_timeidx = try_infer_time_axis(pds)
+
       axes = []
       for i, v in enumerate(VARS):
         full_ds = np.concatenate((hds.variables[v][:,y,x], pds.variables[v][:,y,x]), axis=0)
+        full_idx = h_timeidx.append(p_timeidx).to_pydatetime()
+
         if i > 0:
           ax = plt.subplot(gs[i,0], sharex=axes[0])
           axes.append(ax)
@@ -479,22 +513,25 @@ def climate_ts_plot(args):
       for i, v in enumerate(VARS):
         ax = plt.subplot(gs[i,0])
         with nc.Dataset(os.path.join(args.input_folder, CLIMATE_FILES[0])) as hds:
-          ax.plot(hds.variables[v][:,y, x])
+          hidx = try_infer_time_axis(hds)
+          ax.plot(hidx, hds.variables[v][:,y, x])
           ax.set_title(v)
       plt.suptitle(os.path.join(args.input_folder, CLIMATE_FILES[0]))
       plt.show(block=True)
 
       plt.title(CLIMATE_FILES[1])
       for i, v in enumerate(VARS):
         ax = plt.subplot(gs[i,0])
-        with nc.Dataset(os.path.join(args.input_folder, CLIMATE_FILES[1])) as hds:
-          ax.plot(hds.variables[v][:,y, x])
+        with nc.Dataset(os.path.join(args.input_folder, CLIMATE_FILES[1])) as pds:
+          pidx = try_infer_time_axis(pds)
+          ax.plot(pidx, pds.variables[v][:,y, x])
           ax.set_title(v)
       plt.suptitle(os.path.join(args.input_folder, CLIMATE_FILES[1]))
       plt.show(block=True)
 
   else:
     print("Cmd line arg should be checked such that you can't arrive here.")
+    exit(-1)
 
 
 def tunnel_fast(latvar,lonvar,lat0,lon0):

src/Cohort.cpp

@@ -125,11 +125,16 @@ Cohort::Cohort(int y, int x, ModelData* modeldatapointer):
 Cohort::~Cohort() {
 };
 
+/** Loads data from the projected explicit fire input file into internal datastructures. */
+void Cohort::load_proj_explicit_fire(const std::string& proj_exp_fire_file){
+  fire.load_projected_explicit_data(proj_exp_fire_file, y, x);
+}
+
 /** Provides necessary data to Climate for loading projected climate data*/
 void Cohort::load_proj_climate(const std::string& proj_climate_file){
-
   climate.load_proj_climate(proj_climate_file, y, x);
 }
+
 void Cohort::load_proj_co2(const std::string& proj_co2_file){
   climate.load_proj_co2(proj_co2_file);
 }

src/TEM.cpp

@@ -799,6 +799,7 @@ void advance_model(const int rowidx, const int colidx,
     // Loading projected data instead of historic. FIX?
     runner.cohort.load_proj_climate(modeldata.proj_climate_file);
     runner.cohort.load_proj_co2(modeldata.proj_co2_file);
+    runner.cohort.load_proj_explicit_fire(modeldata.proj_exp_fire_file);
 
     BOOST_LOG_SEV(glg,err) << "MAKE SURE YOUR FIRE INPUTS ARE SETUP CORRECTLY!";
 
@@ -953,7 +954,7 @@ void write_status(const std::string fname, int row, int col, int statusCode) {
   temutil::nc( nc_inq_varid(ncid, "run_status", &statusV) );
 
   // Write data
-  BOOST_LOG_SEV(glg, note) << "WRITING FOR OUTPUT STATUS FOR (row, col): " << row << ", " << col << "\n";
+  BOOST_LOG_SEV(glg, note) << "WRITING OUTPUT STATUS FOR (row, col): " << row << ", " << col << "\n";
   temutil::nc( nc_put_var1_int(ncid, statusV, start, &statusCode) );
 
   /* Close the netcdf file. */

src/TEMUtilityFunctions.cpp

@@ -563,11 +563,16 @@ namespace temutil {
     std::vector<DTYPE> data2;
 
     BOOST_LOG_SEV(glg, debug) << "Grab the data from the netCDF file...";
-    if (the_type == NC_INT64 || NC_INT) {
+    if (the_type == NC_INT) {
       int dataI[timeD_len];
       temutil::nc( nc_get_vara_int(ncid, timeseries_var, start, count, &dataI[0]) );
       unsigned dataArraySize = sizeof(dataI) / sizeof(DTYPE);
       data2.insert(data2.end(), &dataI[0], &dataI[dataArraySize]);
+    } else if (the_type == NC_INT64) {
+      int64_t dataI64[timeD_len];
+      temutil::nc( nc_get_vara(ncid, timeseries_var, start, count, &dataI64[0]) );
+      unsigned dataArraySize = sizeof(dataI64) / sizeof(DTYPE);
+      data2.insert(data2.end(), &dataI64[0], &dataI64[dataArraySize]);
     } else if (the_type == NC_FLOAT) {
       float dataF[timeD_len];
       temutil::nc( nc_get_vara_float(ncid, timeseries_var, start, count, &dataF[0]) );
@@ -1103,6 +1108,8 @@ namespace temutil {
 
   template std::vector<int> get_timeseries<int>(const std::string &filename,
       const std::string &var, const int y, const int x);
+  template std::vector<int64_t> get_timeseries<int64_t>(const std::string &filename,
+      const std::string &var, const int y, const int x);
   template std::vector<float> get_timeseries<float>(const std::string &filename,
       const std::string &var, const int y, const int x);
 

src/WildFire.cpp

@@ -55,7 +55,7 @@ WildFire::WildFire(const std::string& fri_fname,
     this->exp_burn_mask = temutil::get_timeseries<int>(exp_fname, "exp_burn_mask", y, x);
     this->exp_fire_severity = temutil::get_timeseries<int>(exp_fname, "exp_fire_severity", y, x);
     this->exp_jday_of_burn = temutil::get_timeseries<int>(exp_fname, "exp_jday_of_burn", y, x);
-    this->exp_area_of_burn = temutil::get_timeseries<int>(exp_fname, "exp_area_of_burn", y, x);
+    this->exp_area_of_burn = temutil::get_timeseries<int64_t>(exp_fname, "exp_area_of_burn", y, x);
   }//End critical(exp_fir) 
 
   this->slope = cell_slope;
@@ -67,6 +67,14 @@ WildFire::WildFire(const std::string& fri_fname,
 
 }
 
+void WildFire::load_projected_explicit_data(const std::string& exp_fname, int y, int x) {
+    BOOST_LOG_SEV(glg, info) << "Setting up explicit fire data...";
+    this->exp_burn_mask = temutil::get_timeseries<int>(exp_fname, "exp_burn_mask", y, x);
+    this->exp_fire_severity = temutil::get_timeseries<int>(exp_fname, "exp_fire_severity", y, x);
+    this->exp_jday_of_burn = temutil::get_timeseries<int>(exp_fname, "exp_jday_of_burn", y, x);
+    this->exp_area_of_burn = temutil::get_timeseries<int64_t>(exp_fname, "exp_area_of_burn", y, x);
+}
+
 /** Assemble and return a string with a bunch of data from this class */
 std::string WildFire::report_fire_inputs() {