Add postprocessing_example as script and adjust slides

lectures/debugging-strategies/slides.qmd

@@ -542,10 +542,10 @@ issues in line 9 and 13
 ## Postprocessing example (1/4)
 
  - Extract data from netCDF file and calculate interpolation weights
- - File located [here](static/vmro3_input4MIPs_ozone_1850-1855_rev.nc)
+ - Python script located [here](static/postprocessing_example.py), netCDF file [here](static/vmro3_input4MIPs_ozone_1850-1855_rev.nc)
  - Works for the active time stamps, but not for the ones in the comment
 
-```python
+```python{.python startFrom="46"}
 filename = 'vmro3_input4MIPs_ozone_1850-1855_rev.nc'
 timestamp_start = '1851-01-01T00:00:00.000'  # '1850-01-01T00:00:00.000'
 timestamp_end   = '1851-03-31T00:00:00.000'  # '1850-03-31T00:00:00.000'
@@ -607,12 +607,6 @@ def Select_Date(ds, model_date, method):
 
 <!--
 ```python
-def Open_File(filename):
-    import xarray as xr
-
-    return xr.open_dataset(filename)\
-        .convert_calendar('standard', use_cftime=True)
-
 def Select_Date(ds, model_date, method='ffill'):
     import datetime as dt
 

lectures/debugging-strategies/static/postprocessing_example.py

+def Open_File(filename):
+    import xarray as xr
+
+    return xr.open_dataset(filename)\
+        .convert_calendar('standard', use_cftime=True)
+
+def Select_Date(ds, model_date, method):
+    import datetime as dt
+
+    ds_elem = ds.sel(time=model_date, method=method)
+    return dt.datetime.strptime(str(ds_elem['time'].values),
+        '%Y-%m-%d %H:%M:%S')
+
+def Get_Time_Objects(timestamp_start, timestamp_end, timestep_days):
+    import datetime as dt
+
+    time_format = '%Y-%m-%dT%H:%M:%S.%f'
+    start_date = dt.datetime.strptime(timestamp_start, time_format)
+    end_date = dt.datetime.strptime(timestamp_end, time_format)
+    timestep   = dt.timedelta(days=timestep_days)
+    return [start_date, end_date, timestep]
+
+def Calculate_Weights(filename, timestamp_start, timestamp_end, timestep_days):
+    start_date, end_date, timestep = Get_Time_Objects(
+        timestamp_start=timestamp_start, timestamp_end=timestamp_end,
+        timestep_days=10)
+
+    ds = Open_File(filename=filename)
+
+    model_date = start_date
+    while ( model_date < end_date ):
+        o3_prev_date = Select_Date(ds=ds, model_date=model_date, method='ffill')
+        o3_next_date = Select_Date(ds=ds, model_date=model_date, method='bfill')
+
+        if ( o3_next_date == o3_prev_date ):
+            prev_weight = 0.5
+        else :
+            delta_step_sec = (o3_next_date - o3_prev_date).total_seconds()
+            delta_sec = (model_date - o3_prev_date).total_seconds()
+            prev_weight = 1.0 - delta_sec/delta_step_sec
+
+        next_weight  = 1.0 - prev_weight
+        print('weights for', model_date, 'are', prev_weight, 'and', next_weight)
+        model_date = model_date + timestep
+
+filename = 'vmro3_input4MIPs_ozone_1850-1855_rev.nc'
+timestamp_start = '1851-01-01T00:00:00.000'  # '1850-01-01T00:00:00.000'
+timestamp_end   = '1851-03-31T00:00:00.000'  # '1850-03-31T00:00:00.000'
+timestep_days = 10
+
+Calculate_Weights(filename, timestamp_start, timestamp_end, timestep_days)