global cv clim

enstools/feature/identification/african_easterly_waves/all_years.sh

+#!/bin/bash
+
+for i in {1979..2019}
+do
+   for m in {6..10}
+   do
+      python run_identify.py $i $m
+   done
+done
+

enstools/feature/identification/african_easterly_waves/configuration.py

@@ -10,16 +10,17 @@ from datetime import timedelta
 data_lat = (-35, 35)
 data_lon = (-100, 45)
 
-aew_clim_dir = join(expanduser("~") + '/phd/data/aew/clim/')
-cv_data_dir = join(expanduser("~") + '/phd/data/aew/cv/') # reference where ALL cv data is (for clim calc.)
+aew_clim_dir = '/project/meteo/w2w/C3/fischer/belanger/aew_clim/' # join(expanduser("~") + '/phd/data/aew/clim/')
+cv_data_dir = '/project/meteo/w2w/C2/athul/data_athul/AEW_cv_data/' # reference where ALL cv data is (for clim calc.)
+out_dir = '/project/meteo/w2w/C3/fischer/belanger/out/'
+plot_dir = '/project/meteo/w2w/C3/fischer/belanger/plots/'
 
-
-# construct clim file, regenerate for each window?
 def get_clim_file():
-    fn = (aew_clim_dir + "cv_clim_" + str(abs(data_lat[1])) + ('N' if data_lat[1] > 0 else 'S') + "_"
-          + str(abs(data_lat[0])) + ('N' if data_lat[0] > 0 else 'S') + "_"
-          + str(abs(data_lon[0])) + ('E' if data_lon[0] > 0 else 'W') + "_"
-          + str(abs(data_lon[1])) + ('E' if data_lon[1] > 0 else 'W')) + ".nc"
+    fn = aew_clim_dir + "cv_clim_era5.nc"
+    #fn = (aew_clim_dir + "cv_clim_" + str(abs(data_lat[1])) + ('N' if data_lat[1] > 0 else 'S') + "_"
+    #      + str(abs(data_lat[0])) + ('N' if data_lat[0] > 0 else 'S') + "_"
+    #      + str(abs(data_lon[0])) + ('E' if data_lon[0] > 0 else 'W') + "_"
+    #      + str(abs(data_lon[1])) + ('E' if data_lon[1] > 0 else 'W')) + ".nc"
     return fn
 
 
@@ -30,8 +31,9 @@ wave_filter_lon = (-110, 55)
 levels = [70000]  # 700 hPa
 
 # time of interest, if None all
-start_date = None # '2008-08-01T00:00' #  # '2008-08-01T00:00'
-end_date = None # '2008-08-01T00:00' # None # '2008-08-03T00:00'
+# june-oct is AEW season
+start_date = '2008-06-01T00:00' #  # '2008-08-01T00:00'
+end_date = '2008-10-31T00:00' # None # '2008-08-03T00:00'
 
 # Algorithm parameters
 # max u wind (m/s) (0 = only keep west-propagating). Belanger: 2.5; Berry: 0.0

enstools/feature/identification/african_easterly_waves/identification.py

@@ -16,7 +16,7 @@ from skimage.draw import line
 
 class AEWIdentification(IdentificationTechnique):
 
-    def __init__(self, wt_out_file=False, wt_traj_dir=None, cv='cv', **kwargs):
+    def __init__(self, wt_out_file=False, wt_traj_dir=None, cv='cv', year_summer=None, month=None, **kwargs):
         """
         Initialize the AEW Identification.
 
@@ -24,14 +24,27 @@ class AEWIdentification(IdentificationTechnique):
         ----------
         kwargs
         wt_out_file: output the wavetroughs as new and only out-field in 0.5x0.5
+        year_summer: if set, process AEW season (01.06.-31.10.) of given year
+        
         """
 
         import enstools.feature.identification.african_easterly_waves.configuration as cfg
         self.config = cfg  # config
         self.config.out_traj_dir = wt_traj_dir
         self.config.cv_name = cv
-
+        
+        if year_summer is not None:
+            if month is not None:
+                m_str = str(month).zfill(2)
+                self.config.start_date = str(year_summer) + '-' + m_str + '-01T00:00'
+                self.config.end_date = str(year_summer) + '-' + m_str + '-30T00:00'
+            else:
+                self.config.start_date = str(year_summer) + '-06-01T00:00'
+                self.config.end_date = str(year_summer) + '-10-31T00:00'
+        
         self.config.out_wt = wt_out_file
+        if wt_out_file:
+            self.config.sum_over_all = True # TODO maybe more general.
         self.lock_ = threading.Lock()
         pass
 
@@ -51,7 +64,7 @@ class AEWIdentification(IdentificationTechnique):
 
         if os.path.isfile(clim_file):
             cv_clim = xr.open_dataset(clim_file)
-
+            
         else:
             # generate: need all 40y of CV data.
             print("Climatology file not found. Computing climatology...")
@@ -59,7 +72,11 @@ class AEWIdentification(IdentificationTechnique):
             from .climatology import compute_climatology
             cv_clim = compute_climatology(self.config)
             cv_clim.to_netcdf(clim_file)
-
+        
+        cv_clim = cv_clim.sel(
+            **{lat_str: slice(lat_range[0], lat_range[1])},
+            **{lon_str: slice(lon_range[0], lon_range[1])})
+        
         # --------------- SUBSET DATA ACCORDING TO CFG
         start_date_dt = np.datetime64(self.config.start_date) if self.config.start_date is not None else None
         end_date_dt = np.datetime64(self.config.end_date) if self.config.end_date is not None else None
@@ -69,13 +86,10 @@ class AEWIdentification(IdentificationTechnique):
             **{lat_str: slice(lat_range[0], lat_range[1])},
             **{lon_str: slice(lon_range[0], lon_range[1])},
             time=slice(start_date_dt, end_date_dt))
-        if level_str is None:
-            if not 'level' in dataset.coords:
-                print("No level information given in input. Assume 700hPa.")
-                exit(1)
-            dataset = dataset.expand_dims('level')
-            level_str = 'level'
-        else:
+        
+        # dataset = dataset.expand_dims('level')
+        # level_str = 'level'
+        if level_str is not None:
             dataset = dataset.sel(**{level_str: self.config.levels})
 
         # rename cv_clim dimensions to be same as in data.
@@ -238,6 +252,9 @@ class AEWIdentification(IdentificationTechnique):
             level_str = get_vertical_dim(dataset)
             if level_str is not None:
                 dataset = dataset.squeeze(drop=True)
+                
+            if self.config.sum_over_all:
+                dataset['wavetroughs'] = dataset.wavetroughs.sum(dim='time')
 
         # create met3d like trajectories TODO not really working right now...
         if self.config.out_traj_dir:

enstools/feature/identification/african_easterly_waves/plotting.py

@@ -9,6 +9,7 @@ import cartopy.feature as cfeature
 from datetime import datetime
 from enstools.feature.util.data_utils import pb_str_to_datetime
 from pathlib import Path
+import enstools.feature.identification.african_easterly_waves.configuration as cfg
 
 # plots the wave state (all wavetroughs given specific timestep in a set) ts: pb2.Timestep
 def plot_wavetroughs(ts, fig_name, cv=None):
@@ -63,7 +64,7 @@ def plot_timesteps_from_desc(object_desc, cv=None):
     from enstools.feature.util.data_utils import get_subset_by_description
 
     for set_idx, od_set in enumerate(object_desc.sets):
-        fn = "ts_set_" + str(set_idx)
+        fn = cfg.plot_dir + "ts_set_" + str(set_idx)
 
         cv_set = get_subset_by_description(cv, od_set, '2d')
 
@@ -168,8 +169,8 @@ def plot_wt_list(nodes, fn):
     return figure_name
 
 def plot_track_in_ts(track):
-    fn = os.path.join(str(os.path.expanduser('~')),
-                      Path("phd/data/aew/plots/singletrack_"))
+    
+    fn = cfg.plot_dir + 'singletrack_'
 
     per_ts_wts = dict()
     for edge in track.edges:
@@ -184,7 +185,7 @@ def plot_track_in_ts(track):
         plot_wt_list(nodes, fn + time)
 
 
-def plot_tracks_from_graph(graph_desc, ds=None):
+def plot_tracks_from_desc(graph_desc, ds=None):
 
     from enstools.feature.util.data_utils import get_subset_by_description
 
@@ -192,7 +193,7 @@ def plot_tracks_from_graph(graph_desc, ds=None):
 
         # cv_set = get_subset_by_description(ds, od_set, '2d')
         for set_tr, track in enumerate(od_set.tracks):
-            fn = os.path.join(str(os.path.expanduser('~')), Path("phd/data/aew/plots/set_" + str(set_idx) + "_track_" + str(set_tr)))
+            fn = cfg.plot_dir + 'set_' + str(set_idx) + "_track_" + str(set_tr)
             plot_track_from_graph(track, fn, cv=None)
 
     pass

enstools/feature/identification/african_easterly_waves/run_identify.py

@@ -6,24 +6,28 @@ from datetime import timedelta
 from enstools.feature.identification._proto_gen import african_easterly_waves_pb2
 from os.path import expanduser, join
 from enstools.feature.util.graph import DataGraph
-from enstools.feature.identification.african_easterly_waves.plotting import plot_track_in_ts, plot_timesteps_from_desc, plot_tracks_from_graph
+from enstools.feature.identification.african_easterly_waves.plotting import plot_track_in_ts, plot_timesteps_from_desc, plot_tracks_from_desc
+import enstools.feature.identification.african_easterly_waves.configuration as cfg
+import sys
 
 pipeline = FeaturePipeline(african_easterly_waves_pb2, processing_mode='2d')
 
-# data_dir = join(expanduser("~") + '/phd/data/aew/cv/') # 2008cv.nc
-data_dir = join(expanduser("~"), 'phd/data/aew/comp/helene/') # helene_25.nc
-in_file = join(data_dir, 'helene_25.nc') # 2008cv.nc helene_25.nc
-# out_traj_dir = join(data_dir, 'trajectories/')
+data_dir = cfg.cv_data_dir
+in_file = join(data_dir, '*.nc')
+out_dir = cfg.out_dir
 
+if len(sys.argv) > 1:
+    proc_summer_of_year = int(sys.argv[1])
+    if len(sys.argv) > 2:
+        proc_month_of_year = int(sys.argv[2])
+    
 # init AEWIdentification strategy, can take different parameters
-i_strat = AEWIdentification(wt_out_file=False, cv='cv')
+i_strat = AEWIdentification(wt_out_file=True, cv='cv', year_summer=proc_summer_of_year, month=proc_month_of_year)
 t_strat = AEWTracking()
 
 pipeline.set_identification_strategy(i_strat)
 pipeline.set_tracking_strategy(t_strat) # None TODO
 
-# data_path = join(expanduser("~") + '/phd/data/aew/cv/2008cv.nc')
-# data_path = join(expanduser("~") + '/phd/data/aew/comp/helene/helene_25.nc')
 pipeline.set_data_path(in_file)
 
 # execute pipeline
@@ -35,13 +39,24 @@ pipeline.generate_tracks()
 od = pipeline.get_object_desc() # just the identified objects and their connections in object_desc
 
 # out_netcdf_path = data_path + '_streamers.nc'
-out_json_path = in_file[:-3] + 'aew_desc.json'
-out_graph_path = in_file[:-3] + 'aew_graph.json'
-out_dataset_path = in_file[:-5] + '05_wt.nc'
 
-pipeline.save_result(description_type='json', dataset_path=out_dataset_path, description_path=out_json_path) # , description_path=out_json_path, graph_path=out_graph_path
+if len(sys.argv) == 1:
+    out_json_path = out_dir + 'aew_desc.json'
+    out_dataset_path = out_dir + '05_wt.nc'
+elif len(sys.argv) == 2:
+    out_json_path = out_dir + 'aew_desc_' + str(proc_summer_of_year) + '.json'
+    out_dataset_path = out_dir + '05_wt_' + str(proc_summer_of_year) + '.nc'
+else:
+    m_str = str(proc_month_of_year).zfill(2)
+    out_json_path = out_dir + 'aew_desc_' + str(proc_summer_of_year) + '_' + m_str + '.json'
+    out_dataset_path = out_dir + '05_wt_' + str(proc_summer_of_year) + '_' + m_str + '.nc'
+    
+
+
+pipeline.save_result(description_type='json', description_path=out_json_path, dataset_path=out_dataset_path) # dataset_path=out_dataset_path, 
+# , description_path=out_json_path, graph_path=out_graph_path
+
+# print("Plot.")
+# plot_timesteps_from_desc(od, pipeline.get_data())
+# plot_tracks_from_desc(od, None)
 
-print("Plot.")
-#plot_timesteps_from_desc(od, pipeline.get_data())
-# plot_tracks_from_graph(od, None)
-# plot_track_in_ts(od.sets[0].tracks[8])
\ No newline at end of file