diff --git a/enstools/feature/identification/african_easterly_waves/configuration.py b/enstools/feature/identification/african_easterly_waves/configuration.py index 25356238944f8cd4505aedec4d3e2bfe7b1f50ce..a1856f19e18aedc60d54aeef9933ffe821713431 100644 --- a/enstools/feature/identification/african_easterly_waves/configuration.py +++ b/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') diff --git a/enstools/feature/identification/african_easterly_waves/identification.py b/enstools/feature/identification/african_easterly_waves/identification.py index 9b940624e0c454d3e4de61f8ccf774e7499417f0..ccfef8a9e080dc4a1f690011f544e6ed06284c2b 100644 --- a/enstools/feature/identification/african_easterly_waves/identification.py +++ b/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])}, diff --git a/enstools/feature/identification/african_easterly_waves/plotting.py b/enstools/feature/identification/african_easterly_waves/plotting.py index bdc22294194ef5afa97b4d01d5ab279ac105195f..a69843c14f83a48dd75f358d3d097145abf825cb 100644 --- a/enstools/feature/identification/african_easterly_waves/plotting.py +++ b/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) diff --git a/enstools/feature/identification/african_easterly_waves/run_identify.py b/enstools/feature/identification/african_easterly_waves/run_identify.py index 82b881bd37a6acd5f8b00b01fb08a571ced24227..0da89df5053ddeb76dbf0417a5f766cb1883817c 100644 --- a/enstools/feature/identification/african_easterly_waves/run_identify.py +++ b/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'