icon support, moved related dirs

enstools/feature/identification/african_easterly_waves/configuration.py

@@ -11,13 +11,24 @@ data_lat = (-35, 35)
 data_lon = (-100, 45)
 
 aew_clim_dir = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/era5/cv_clim_era5.nc' # TODO # '/home/christoph/phd/data/framework_example_ds/aew/' # '/project/meteo/w2w/C3/fischer/belanger/aew_clim/' # 
+
 aew_kitweather_ecmwf_dir = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/ecmwf-hres/forecasts/'
-latest_run_info_file = '/home/iconeps/plots/latest_aew_run.txt'
+num_files_ecmwf = 41 # files in finished forecast data
+file_prefix_ecmwf = 'ecmwf-hres_latlon_1.0deg_'
+file_prefix_ecmwf_rain = 'ecmwf-hres_latlon_0.4deg_'
+plot_dir_ecmwf = '/home/iconeps/plots/aew/ecmwf-hres/' # '/project/meteo/w2w/C3/fischer/belanger/plots/'
+
+aew_kitweather_icon_dir = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/icon-global-det/forecasts/'
+num_files_icon = 31 # files in finished forecast data
+file_prefix_icon = 'icon-global-det_latlon_1.0deg_'
+file_prefix_icon_rain = 'icon-global-det_latlon_0.25deg_'
+plot_dir_icon = '/home/iconeps/plots/aew/icon-global-det/' # '/project/meteo/w2w/C3/fischer/belanger/plots/'
+
+latest_run_info_file_name = 'latest_run.txt'
 
 in_files = '/home/christoph/phd/data/aew/cv/2008cv.nc' # '/home/christoph/phd/data/framework_example_ds/aew/cv_aug_08.nc'
 
 # out_dir = join(expanduser("~") + '/phd/data/aew/out/') # '/project/meteo/w2w/C3/fischer/belanger/out/'
-plot_dir = join('/home/iconeps/plots/aew/') # '/project/meteo/w2w/C3/fischer/belanger/plots/'
 
 cv_data_dir = '/project/meteo/w2w/C2/athul/data_athul/AEW_cv_data/' # reference where ALL cv data is (only for clim calc.)
 cv_data_ex = join(cv_data_dir, '*.nc')

enstools/feature/identification/african_easterly_waves/identification.py

@@ -80,8 +80,12 @@ class AEWIdentification(IdentificationTechnique):
         # --------------- 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
-        # get the data we want to investigate
+        
+        # if data is lon=0..360, change it to -180..180
+        dataset.coords[lon_str] = (dataset.coords[lon_str] + 180) % 360 - 180
+        dataset = dataset.sortby(dataset[lon_str])
 
+        # get the data we want to investigate
         dataset = dataset.sel(
             **{lat_str: slice(lat_range[0], lat_range[1])},
             **{lon_str: slice(lon_range[0], lon_range[1])},

enstools/feature/identification/african_easterly_waves/plotting.py

@@ -303,10 +303,10 @@ def plot_differences(set_graph, tracks, ds=None, tp=None, plot_prefix=None):
     if plot_prefix is None:
         plot_prefix = cfg.plot_dir
     # create subdirs if needed
-    os.makedirs(plot_prefix, exist_ok=True)
+    # os.makedirs(plot_prefix, exist_ok=True)
 
     for date in dates_list:
-        fig_name = plot_prefix + "aew_" + date[0:4] + date[5:7] + date[8:10] + "T" + date[11:13] + ".png"
+        fig_name = plot_prefix + date[0:4] + date[5:7] + date[8:10] + "T" + date[11:13] + ".png"
         ds_ss = ds.sel(time=date)
         try:
             tp_ss = tp.sel(time=date)

enstools/feature/identification/african_easterly_waves/run_identify.py

@@ -35,9 +35,22 @@ if len(sys.argv) == 3 and sys.argv[1] == '-kw' and sys.argv[2] == 'ana':
     in_file = data_fc_root + "*/*000h_tropicalvars.nc"
     print("Executing for: " + in_file)
 
-elif len(sys.argv) >= 3 and sys.argv[1] == '-kw' and sys.argv[2] == 'ecmwf_fc':
+elif len(sys.argv) >= 3 and sys.argv[1] == '-kw':
+    print("Executing in kitweather mode...")
     # kitweather: use last 7 days of analysis and the ecmwf forecast
-    data_fc_root = cfg.aew_kitweather_ecmwf_dir
+    if sys.argv[2] == 'ecmwf_fc':
+        data_fc_root = cfg.aew_kitweather_ecmwf_dir
+        num_files_in_fc = cfg.num_files_ecmwf
+        file_prefix = cfg.file_prefix_ecmwf
+        file_prefix_rain = cfg.file_prefix_ecmwf_rain
+    elif sys.argv[2] == 'icon_fc':
+        data_fc_root = cfg.aew_kitweather_icon_dir
+        num_files_in_fc = cfg.num_files_icon
+        file_prefix = cfg.file_prefix_icon
+        file_prefix_rain = cfg.file_prefix_icon_rain
+    else:
+        print("Unknown command line parameter " + sys.argv[2] + ", expected ecmwf_fc or icon_fc.")
+        exit(1)
 
     # get latest subdirectory time as what we choose
     all_subdirs = glob.glob(data_fc_root + "*/")
@@ -45,50 +58,47 @@ elif len(sys.argv) >= 3 and sys.argv[1] == '-kw' and sys.argv[2] == 'ecmwf_fc':
     # if as additional argument run_YYYYMMDDHH, use this, otherwise find newest.
     if len(sys.argv) == 4:
         latest_subdir = data_fc_root + sys.argv[3] + '/'
-        forecast_files_glob = latest_subdir + "ecmwf-hres_latlon_1.0deg_*"
+        forecast_files_glob = latest_subdir + file_prefix + "*"
         fc_file_list = glob.glob(forecast_files_glob)
 
-        fc_rain_file_list = [fc_file.replace("1.0", "0.4").replace("tropicalvars", "tp") for fc_file in fc_file_list if "000h" not in fc_file]
-        if len(fc_file_list) < 41:
+        fc_rain_file_list = [fc_file.replace(file_prefix, file_prefix_rain).replace("tropicalvars", "tp") for fc_file in fc_file_list if "000h" not in fc_file]
+        if len(fc_file_list) < num_files_in_fc:
             print("Missing files. Found " + str(len(fc_file_list)) + " files in forecast.")
             print("Exit.")
             exit(1)
     else:
-    
-        print(data_fc_root + "*/")
 
         print("Collecting forecast files...")
         sorted_subdirs = sorted(all_subdirs)
         latest_subdir = sorted_subdirs[-1]
-        forecast_files_glob = latest_subdir + "ecmwf-hres_latlon_1.0deg_*"
+        forecast_files_glob = latest_subdir + file_prefix + "*"
         fc_file_list = glob.glob(forecast_files_glob)
 
-        fc_rain_file_list = [fc_file.replace("1.0", "0.4").replace("tropicalvars", "tp") for fc_file in fc_file_list if "000h" not in fc_file]
+        fc_rain_file_list = [fc_file.replace(file_prefix, file_prefix_rain).replace("tropicalvars", "tp") for fc_file in fc_file_list if "000h" not in fc_file]
 
-        if len(fc_file_list) < 41:
+        if len(fc_file_list) < num_files_in_fc:
             print("Expected 41 files in " + forecast_files_glob + ", got " + str(len(fc_file_list)))
             print("Trying previous timestep...")
             latest_subdir = sorted_subdirs[-2]
-            forecast_files_glob = latest_subdir + "ecmwf-hres_latlon_1.0deg_*"
+            forecast_files_glob = latest_subdir + file_prefix + "*"
             fc_file_list = glob.glob(forecast_files_glob)
-            fc_rain_file_list = [fc_file.replace("1.0", "0.4").replace("tropicalvars", "tp") for fc_file in fc_file_list if "000h" not in fc_file]
+            fc_rain_file_list = [fc_file.replace(file_prefix, file_prefix_rain).replace("tropicalvars", "tp") for fc_file in fc_file_list if "000h" not in fc_file]
 
-            if len(fc_file_list) < 41:
+            if len(fc_file_list) < num_files_in_fc:
                 print("Missing files as well. Found " + str(len(fc_file_list)) + " files.")
                 print("Exit.")
                 exit(1)
-    print("Found all 41 forecast files at " + forecast_files_glob + ".")
+    print("Found all " + str(num_files_in_fc) + " forecast files at " + forecast_files_glob + ".")
 
     # get last 7 days of analysis: 000h from previous runs
     
     print("Collecting analysis files...")
-    data_fc_root = cfg.aew_kitweather_ecmwf_dir
     all_subdirs_by_time = sorted(all_subdirs)
 
     # get analysis files starting at latest_subdirs 7 days in the past
     latest_subdir_rel_index = all_subdirs_by_time.index(latest_subdir) - len(all_subdirs_by_time)
     last_7d_ana_subdirs = all_subdirs_by_time[latest_subdir_rel_index-28:latest_subdir_rel_index] # last 28 timesteps = last 7 days
-    last_7d_ana_glob = [sd + "ecmwf-hres_latlon_1.0deg_*_000h_tropicalvars.nc" for sd in last_7d_ana_subdirs]
+    last_7d_ana_glob = [sd + file_prefix + "*_000h_tropicalvars.nc" for sd in last_7d_ana_subdirs]
 
     ana_file_list = []
     for ana_ts_glob in last_7d_ana_glob:
@@ -104,7 +114,8 @@ elif len(sys.argv) >= 3 and sys.argv[1] == '-kw' and sys.argv[2] == 'ecmwf_fc':
 
     print("Collecting rain data...")
     # just get all tp 0.4deg 6h
-    ana_rain_files = glob.glob(data_fc_root + "*/ecmwf-hres_latlon_0.4deg_*_006h_tp.nc")
+    ana_rain_files = glob.glob(data_fc_root + "*/" + file_prefix_rain + "*_006h_tp.nc")
+
 
     # load forecast files separately: need to compute deltas from tp.
     rain_fc_ds = xr.open_mfdataset(fc_rain_file_list)
@@ -112,6 +123,7 @@ elif len(sys.argv) >= 3 and sys.argv[1] == '-kw' and sys.argv[2] == 'ecmwf_fc':
     # rain_fc_tp_diff = rain_fc_tp.differentiate(coord="time", datetime_unit="6h")
 
     times = rain_fc_tp.time.values
+
     for t_idx, time in reversed(list(enumerate(times))):
         if t_idx > 0:
             rain_fc_tp.loc[dict(time=time)] = rain_fc_tp.isel(time=t_idx) - rain_fc_tp.isel(time=(t_idx-1))
@@ -129,6 +141,7 @@ elif len(sys.argv) >= 3 and sys.argv[1] == '-kw' and sys.argv[2] == 'ecmwf_fc':
 
     print(rain_ana_ds.time.data)
     print(rain_fc_ds.time.data)
+
     # select analysis data up to forecast start. sometimes overlap if future analysis gets downloaded. done above.
 
     all_rain_ds = xr.merge([rain_ana_ds, rain_fc_ds])
@@ -136,7 +149,12 @@ elif len(sys.argv) >= 3 and sys.argv[1] == '-kw' and sys.argv[2] == 'ecmwf_fc':
     # change total precipitation to hourly precip
     all_rain_ds.tp.attrs['units'] = 'mm hr-1'
     all_rain_ds.tp.attrs['long_name'] = 'Precipitation rate'
-    all_rain_ds['tp'] = all_rain_ds.tp / 6.0 * 1000.0 # from 6hrly (downloaded) to hourly rate and m to mm
+    # kg m-2 = mm /// 1000 mm = m
+    # ecmwf is m, icon is kg m-2
+    # from 6hrly (downloaded) to hourly rate and m to mm
+    all_rain_ds['tp'] = all_rain_ds.tp / 6.0
+    if sys.argv[2] == 'ecmwf_fc':
+        all_rain_ds['tp'] = all_rain_ds.tp / 6.0 * 1000.0
 
     print("Done collecting files.")
 
@@ -177,9 +195,20 @@ for trackable_set in od.sets:
     ds = pipeline.get_data()
     ds_set = get_subset_by_description(ds, trackable_set, '2d')
 
-    if len(sys.argv) >= 3 and sys.argv[1] == '-kw' and sys.argv[2] == 'ecmwf_fc':
+    if len(sys.argv) >= 3 and sys.argv[1] == '-kw':
         time_dir = os.path.basename(os.path.normpath(latest_subdir))
-        plot_differences(g, tracks, ds=ds_set, tp=all_rain_ds, plot_prefix=cfg.plot_dir + time_dir + "/")
+        if sys.argv[2] == 'ecmwf_fc':
+            plot_dir = cfg.plot_dir_ecmwf
+            plot_prefix = cfg.plot_dir_ecmwf + time_dir + "/aew_ecmwf-hres_"
+        elif sys.argv[2] == 'icon_fc':
+            plot_dir = cfg.plot_dir_icon
+            plot_prefix = cfg.plot_dir_icon + time_dir + "/aew_icon-global-det_"
+        else:
+            print("Unknown " + sys.argv[2])
+            exit(1)
+        plot_differences(g, tracks, ds=ds_set, tp=all_rain_ds, plot_prefix=plot_prefix)
+
+
     elif len(sys.argv) >= 3 and sys.argv[1] == '-kw' and sys.argv[2] == 'ana':
         plot_differences(g, tracks, ds=ds_set, tp=rain_ds, plot_prefix=cfg.plot_dir + "ana/")
     else:
@@ -189,28 +218,30 @@ time_dir = os.path.basename(os.path.normpath(latest_subdir))
 # no out data besides plots on kitweather
 if sys.argv[1] == '-kw':
     # delete old plots
-    subdirs = [dI for dI in os.listdir(cfg.plot_dir) if os.path.isdir(os.path.join(cfg.plot_dir,dI))]
+    subdirs = [dI for dI in os.listdir(plot_dir) if os.path.isdir(os.path.join(plot_dir,dI))]
     for sd in subdirs: # for each subdir in plot dir
 
-        if not sd == time_dir and datetime.fromtimestamp(os.path.getmtime(os.path.join(cfg.plot_dir, sd))) < datetime.now() - timedelta(days=7): # not touched in a week? delete it.
-            print("Removing directory " + str(os.path.join(cfg.plot_dir, sd)))
-            shutil.rmtree(os.path.join(cfg.plot_dir, sd))
+        if not sd == time_dir and datetime.fromtimestamp(os.path.getmtime(os.path.join(plot_dir, sd))) < datetime.now() - timedelta(days=7): # not touched in a week? delete it.
+            print("Removing directory " + str(os.path.join(plot_dir, sd)))
+            shutil.rmtree(os.path.join(plot_dir, sd))
 
 
     # All done. Update text file containing time of latest finished run.
     yyyymmddhh = time_dir[4:]
 
-    with open(cfg.latest_run_info_file, 'w+') as info_file:
+    with open(plot_dir + cfg.latest_run_info_file_name, 'w+') as info_file:
         info_file.write(yyyymmddhh)
         
     # All done. scp data over to webserver.
     # TODO when port 22 free test.
     # TODO improve plotting performance
+    # TODO depending on ecmwf or icon
+    """
     path_webserver = '/home/iconeps/Data3/plots/ecmwf/aew_prediction_maps/'
     print('scp -r ' + cfg.plot_dir + time_dir + '/ ' + 'iconeps@imk-tss-web.imk-tro.kit.edu:' + path_webserver)
     os.system('scp -r ' + cfg.plot_dir + time_dir + '/ ' + 'iconeps@imk-tss-web.imk-tro.kit.edu:' + path_webserver)
     os.system('scp ' + cfg.latest_run_info_file + ' iconeps@imk-tss-web.imk-tro.kit.edu:' + path_webserver)
-
+    """
     exit()
 
 # out_netcdf_path = data_path + '_streamers.nc'