Update spinup args

Author: btobers

pygem/bin/run/run_inversion.py

@@ -280,13 +280,15 @@ def main():
         nargs='+',
         help='Randoph Glacier Inventory glacier number (can take multiple)',
     )
-    parser.add_argument(
-        '-rgi_glac_number_fn',
-        action='store',
-        type=str,
-        default=None,
-        help='Filepath containing list of rgi_glac_number, helpful for running batches on spc',
-    ),
+    (
+        parser.add_argument(
+            '-rgi_glac_number_fn',
+            action='store',
+            type=str,
+            default=None,
+            help='Filepath containing list of rgi_glac_number, helpful for running batches on spc',
+        ),
+    )
     parser.add_argument(
         '-calibrate_regional_glen_a',
         type=str2bool,
@@ -298,13 +300,13 @@ def main():
         type=float,
         default=None,
         help="User-selected inversion Glen's creep parameter value",
-    )  
+    )
     parser.add_argument(
         '-fs',
         type=float,
         default=None,
         help="User-selected inversion Orleam's sliding factor value",
-    )    
+    )
     parser.add_argument(
         '-ncores',
         action='store',
@@ -354,7 +356,7 @@ def main():
         run,
         ncores=args.ncores,
         calibrate_regional_glen_a=args.calibrate_regional_glen_a,
-        glen_a=args.glen_a, 
+        glen_a=args.glen_a,
         fs=args.fs,
         reset_gdirs=args.reset_gdirs,
         debug=args.debug,

pygem/bin/run/run_spinup.py

@@ -30,10 +30,10 @@
 def run(glacno_list, mb_model_params, debug=False, **kwargs):
     # remove None kwargs
     kwargs = {k: v for k, v in kwargs.items() if v is not None}
-
     main_glac_rgi = modelsetup.selectglaciersrgitable(glac_no=glacno_list)
     # model dates
-    dt = modelsetup.datesmodelrun(startyear=1979, endyear=2019)
+    dt = modelsetup.datesmodelrun(startyear=2000 - kwargs.get('spinup_period', 20), endyear=kwargs.get('ye', 2020))
+
     # load climate data
     ref_clim = class_climate.GCM(name='ERA5')
 
@@ -186,7 +186,7 @@ def main():
         help='Fixed spinup period (years). If not provided, OGGM default is used.',
     )
     parser.add_argument('-target_yr', type=int, default=None)
-    parser.add_argument('-ye', type=int, default=2020)
+    parser.add_argument('-ye', type=int, default=None)
     parser.add_argument(
         '-ncores',
         action='store',
@@ -242,6 +242,7 @@ def main():
         mb_model_params=args.mb_model_params,
         target_yr=args.target_yr,
         spinup_period=args.spinup_period,
+        ye=args.ye,
     )
     # parallel processing
     print(f'Processing with {ncores} cores... \n{glac_no_lsts}')

pygem/mcmc.py

@@ -230,15 +230,15 @@ def log_likelihood(self, m):
             if torch.all(pred == float('-inf')):
                 # Invalid model output -> assign -inf likelihood
                 return torch.tensor([-float('inf')])
-        
+
             if i == 0:
                 # --- Base case: mass balance likelihood ---
                 log_likehood += log_normal_density(
                     self.obs[i][0],  # observed values
                     mu=pred,  # predicted values
                     sigma=self.obs[i][1],  # observation uncertainty
                 )
-            
+
             elif i == 1 and len(m) > 3:
                 # --- Extended case: apply density scaling to get binned elevation change ---
                 # Create density field, separate values for ablation/accumulation zones
@@ -251,9 +251,9 @@ def log_likelihood(self, m):
                 )  # scale prediction by model density values (convert from m ice to m thickness change)
 
                 log_likehood += log_normal_density(
-                    self.obs[i][0],         # observations
-                    mu=pred,                # scaled predictions
-                    sigma=self.obs[i][1],   # uncertainty
+                    self.obs[i][0],  # observations
+                    mu=pred,  # scaled predictions
+                    sigma=self.obs[i][1],  # uncertainty
                 )
         return log_likehood
 

pygem/plot/graphics.py

@@ -296,16 +296,17 @@ def plot_resid_histogram(obs, preds, title, fontsize=8, show=False, fpath=None):
     plt.close(fig)
 
 
-def plot_mcmc_elev_change_1d(preds, fls, obs, ela, title, fontsize=8, rate=True, uniform_area=True, show=False, fpath=None):
-
+def plot_mcmc_elev_change_1d(
+    preds, fls, obs, ela, title, fontsize=8, rate=True, uniform_area=True, show=False, fpath=None
+):
     bin_z = np.array(obs['bin_centers'])
     bin_edges = np.array(obs['bin_edges'])
 
     # get initial thickness and surface area
     initial_area = fls[0].widths_m * fls[0].dx_meter
     initial_thickness = getattr(fls[0], 'thick', None)
     initial_surface_h = getattr(fls[0], 'surface_h', None)
-    # sort initial surface height 
+    # sort initial surface height
     sorting = np.argsort(initial_surface_h)
     initial_surface_h = initial_surface_h[sorting]
     initial_area = initial_area[sorting]
@@ -319,7 +320,7 @@ def plot_mcmc_elev_change_1d(preds, fls, obs, ela, title, fontsize=8, rate=True,
         initial_area = obs['bin_area']
 
     if uniform_area:
-        xvals = np.nancumsum(initial_area)*1e-6
+        xvals = np.nancumsum(initial_area) * 1e-6
     else:
         xvals = bin_z
 
@@ -365,7 +366,6 @@ def elev_to_cum_area(x):
         ax[t].axhline(y=0, c='grey', lw=0.5)
         preds = np.stack(preds)
 
-
         ax[t].fill_between(
             xvals,
             (obs['dh'][:, t] - obs['dh_sigma'][:, t]) / nyrs[t],

pygem/shop/elevchange1d.py

@@ -80,7 +80,7 @@ def elev_change_1d_to_gdir(
         - Each element in 'dates' defines one elevation change period with a start and end date,
         stored as strings in 'YYYY-MM-DD' format.
         - 'bin_edges' should be a list length N containing the elevation values of each bin edge.
-        - 'bin_area' should be a list of length N-1 containing the bin areas given by the 'ref_dem' (optional).        
+        - 'bin_area' should be a list of length N-1 containing the bin areas given by the 'ref_dem' (optional).
         - 'dh' should contain M lists of length N-1, where M is the number of periods and N is the number of bin edges. Units are in meters.
         - 'dh_sigma' should either be M lists of length N-1 (matching 'dh') or a single scalar value. Units are in meters.
         - Each list in 'dh' (and optionally 'dh_sigma') corresponds exactly to one period in 'dates'.
@@ -216,14 +216,14 @@ def validate_elev_change_1d_structure(data):
     ref_dem_year = data['ref_dem_year']
     if not isinstance(ref_dem_year, int):
         raise ValueError("'ref_dem_year' must be an integer year.")
-    
+
     # Validate bin_area if present
     if 'bin_area' in data:
         bin_area = data['bin_area']
         if len(bin_area) != N:
             raise ValueError(f"'bin_area' must be a list of length {N}.")
         if not all(isinstance(x, (int, float)) for x in bin_area):
-            raise ValueError(f"All 'bin_area' values must be numeric.")
+            raise ValueError("All 'bin_area' values must be numeric.")
 
     return True
 

pygem/shop/loso25icebridge.py

@@ -515,7 +515,7 @@ def save_elevchange1d(
         note: 'ref_dem' is the reference DEM used for elevation-binning.
         'ref_dem_year' is the acquisition year of the reference DEM.
         'dates' are tuples (or length-2 sublists) of the start and stop date of an individual elevation change record
-        and are stored as strings in 'YYYY-MM-DD' format. 
+        and are stored as strings in 'YYYY-MM-DD' format.
         'bin_edges' should be a list length N containing the elevation values of each bin edge.
         'bin_area' should be a list of length N-1 containing the bin areas given by the 'ref_dem' (optional).
         'dh' should be M lists of length N-1, where M is the number of date pairs and N is the number of bin edges.