Allow user to read in exogenous scenario in demeter.

demeter/config_reader.py

@@ -175,6 +175,13 @@ def __init__(self, params):
         self.use_constraints = self.valid_integer(run_params.get('use_constraints', 1), 'use_constraints', [0, 1])
         self.spatial_resolution = self.valid_limit(run_params.get('spatial_resolution', 0.25), 'spatial_resolution', [0.0, 1000000.0], 'float')
         self.regrid_resolution = self.valid_limit(run_params.get('regrid_resolution', self.spatial_resolution), 'regrid_resolution', [0.0, 1000000.0], 'float')
+
+        self.external_scenario_PFT_name= (run_params.get('PFT_name_to_replace','Urban'))
+        self.external_scenario = (run_params.get('ext_scenario','SSP1'))
+        self.stitch_external = self.valid_integer(run_params.get('stitch_external', 0), 'stitch_external', [1, 0])
+        self.path_to_external = (run_params.get('path_to_external',self.input_dir))
+
+
         self.errortol = self.valid_limit(run_params.get('errortol', 0.001), 'errortol', [0.0, 1000000.0], 'float')
         self.timestep = self.valid_limit(run_params.get('timestep', 5), 'timestep', [1, 1000000], 'int')
         self.proj_factor = self.valid_limit(run_params.get('proj_factor', 1000), 'proj_factor', [1, 10000000000], 'int')

demeter/demeter_io/writer.py

@@ -12,6 +12,7 @@
 from scipy import io as sio
 import pandas as pd
 from demeter import ncdf_conversion as nc
+import xarray as xr
 
 
 def array_to_csv(arr, out_file):
@@ -35,7 +36,8 @@ def save_array(arr, out_file):
 
 
 def lc_timestep_csv(c, yr, final_landclasses, spat_coords, metric_id_array, gcam_regionnumber, spat_water, cellarea,
-                    spat_ludataharm, metric, units='fraction', write_outputs=False, write_ncdf=False,sce="default",resolution=0.05,write_csv=False,regrid_res=0.05):
+                    spat_ludataharm, metric, units='fraction', write_outputs=False, write_ncdf=False,sce="default",resolution=0.05,write_csv=False,regrid_res=0.05,
+                    stitch_external=0,path_to_external="",external_scenario_PFT_name="",external_scenario=""):
     """Save land cover data for each time step as a CSV file."""
 
     # create out path and file name
@@ -79,14 +81,78 @@ def lc_timestep_csv(c, yr, final_landclasses, spat_coords, metric_id_array, gcam
         write_ncdf = False
 
     if write_ncdf:
+        t_joined= pd.DataFrame(data=arr, columns=columns)
+        if stitch_external==1:
+
+           path_to_sce= path_to_external+external_scenario+str(yr)+".nc"
+
+           if os.path.isfile(path_to_sce):
+               print("User has elected to stitch new scenario with current demeter scenarios. Please make sure all scenarios are stored in correct location")
+               src_raster = xr.open_dataset(path_to_sce)
+               target_resolution = (resolution, resolution)
+               resampled = src_raster
+               longitude_list = generate_scaled_coordinates_helper(-180, 180,resolution, ascending=True)
+
+               # positive to negative scaling required for latitude output
+               latitude_list = generate_scaled_coordinates_helper(-90, 90, resolution, ascending=False)
+
+
+               resampled = resampled.interp(lat=latitude_list,
+                               lon=longitude_list)
+
+
+               sce_df = resampled.to_dataframe()
+
+               # Reset the index to convert the multi-index into columns
+               sce_df.reset_index(inplace=True)
+               # Resample the raster to the target resolution
+               sce_df.dropna(inplace=True)
+               sce_df = sce_df[['lat', 'lon', 'value']]
+
+
+
+               sce_df = sce_df.rename(columns={'value': 'frac','lat':'latitude',
+                                               'lon':'longitude'})
+
+               sce_df.latitude = sce_df.latitude.round(3)
+               sce_df.longitude = sce_df.longitude.round(3)
+
+               #sce_df.to_csv('C:/Projects/sce.csv', index=False)
+               or_df = pd.DataFrame(data=arr, columns=columns)
+
+               or_df.latitude = or_df.latitude.round(3)
+               or_df.longitude = or_df.longitude.round(3)
+
+               #or_df.to_csv('C:/Projects/original.csv', index=False)
+               t_temp = pd.merge(or_df, sce_df, on=['latitude', 'longitude'], how='left')
+
+               t_temp['frac'] = np.where(t_temp['frac'].isna(), 0, t_temp['frac'])
+               t_temp['frac'] = np.where(t_temp['water'] >=0.25, 0, t_temp['frac'])
+               t_temp['adj'] = 1 - t_temp['frac']
+
+               t_temp.dropna(inplace=True)
+
+               lu_file_for_col = or_df.drop(['latitude', 'longitude',external_scenario_PFT_name,'region_id','basin_id'], axis=1)
+               pft_columns = lu_file_for_col.columns
+               t_temp['sum_non_ext']= t_temp[pft_columns].sum(axis=1)
+               t_temp[pft_columns] = t_temp[pft_columns].multiply(t_temp['adj'], axis="index")
+               t_temp[pft_columns] = t_temp[pft_columns].divide(t_temp['sum_non_ext'], axis="index")
+
+               t_temp[external_scenario_PFT_name] = t_temp['frac']
+               t_joined = t_temp.drop(['frac', 'adj','sum_non_ext'], axis=1)
+               t_joined = t_joined.dropna()
+               t_joined = t_joined.reset_index(drop=True)
+               #t_joined.to_csv("C:/Projects/sce.csv")
+
+
         x= nc.DemeterToNetcdf(scenario_name= str(sce),
                        project_name="",
                        start_year=2005,
                        end_year=2005,
                        resolution= resolution,
                        csv_input=write_csv,
                        regrid_resolution=regrid_res,
-                       df=pd.DataFrame(data=arr, columns=columns))
+                       df=t_joined)
 
         x.process_output(input_file_directory=c.lu_csv_output_dir,
                          output_file_directory=c.lu_netcdf_output_dir,
@@ -386,3 +452,35 @@ def max_ascii_rast(arr, out_dir, step, alg='max', nodata=-9999, xll=-180, yll=-9
 
     # create output raster
     arr_to_ascii(final_arr, out_rast, xll=-xll, yll=-yll, cellsize=cellsize, nodata=nodata)
+
+def generate_scaled_coordinates_helper(coord_min: float,
+                                    coord_max: float,
+                                    res: float,
+                                    ascending: bool = True,
+                                    decimals: int = 3) -> np.array:
+        """Generate a list of evenly-spaced coordinate pairs for the output grid based on lat, lon values.
+
+        :param coord_min:                   Minimum coordinate in range.
+        :type coord_min:                    float
+
+        :param coord_max:                   Maximum coordinate in range.
+        :type coord_max:                    float
+
+        :param ascending:                   Ascend coordinate values if True; descend if False
+        :type ascending:                    bool
+
+        :param decimals:                    Number of desired decimals to round to.
+        :type decimals:                     int
+
+        :returns:                           Array of coordinate values.
+
+        """
+
+        # distance between centroid and edge of grid cell
+        center_spacing = res / 2
+
+        if ascending:
+            return np.arange(coord_min + center_spacing, coord_max, res).round(decimals)
+
+        else:
+            return np.arange(coord_max - center_spacing, coord_min, -res).round(decimals)

demeter/process.py

@@ -37,6 +37,10 @@ def __init__(self, config, s, step_idx, step, write_outputs=False):
         self.sce = self.config.scenario
         self.res = self.config.spatial_resolution
         self.regrid_res = self.config.regrid_resolution
+        self.stitch_external= self.config.stitch_external
+        self.path_to_external= self.config.path_to_external
+        self.external_scenario_PFT_name = self.config.external_scenario_PFT_name
+        self.external_scenario= self.config.external_scenario
         # populate
         self.output_df = self.process()
 
@@ -145,7 +149,8 @@ def outputs(self):
                 write_ncdf =True
             return wdr.lc_timestep_csv(self.config, self.step, self.s.final_landclasses, self.s.spat_coords, orig_spat_aez,
                                 self.s.spat_region, self.s.spat_water, self.s.cellarea, self.s.spat_ludataharm,
-                                self.config.metric, self.config.tabular_units, self.write_outputs, write_ncdf, self.sce, self.res, write_csv, self.regrid_res)
+                                self.config.metric, self.config.tabular_units, self.write_outputs, write_ncdf, self.sce, self.res, write_csv, self.regrid_res,
+                                       self.stitch_external,self.path_to_external,self.external_scenario_PFT_name,self.external_scenario)
 
     def process(self):
         """