diff --git a/enstools/feature/identification/african_easterly_waves/climatology.py b/enstools/feature/identification/african_easterly_waves/climatology.py index 8bd6ed484a5c61b2b485e006d38a928e52734b3d..964b83ec97ca1a35568916186cd2044ce67d4c87 100644 --- a/enstools/feature/identification/african_easterly_waves/climatology.py +++ b/enstools/feature/identification/african_easterly_waves/climatology.py @@ -26,6 +26,7 @@ def compute_climatology(cfg): last = '2019-12-31' cv_database_dir = cfg.cv_data_dir + print(cv_database_dir) f = xr.open_mfdataset(cv_database_dir + "*.nc").sel( # plev=slice(plev, plev), diff --git a/enstools/feature/identification/african_easterly_waves/configuration.py b/enstools/feature/identification/african_easterly_waves/configuration.py index 2a478d122bb5c306def9b3b68ebd3378dc8f4100..62388c78f96bdb703f7247a5ab7bc4116c7e017f 100644 --- a/enstools/feature/identification/african_easterly_waves/configuration.py +++ b/enstools/feature/identification/african_easterly_waves/configuration.py @@ -11,25 +11,12 @@ import numpy as np data_lat = (-15, 35) data_lon = (-100, 45) -aew_clim_dir = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/era5/cv_clim_era5.nc' # 'C:\\Users\\Christoph\\phd\\data\\enstools-feature\\cv_clim_era5.nc' # '/home/christoph/phd/data/aew/clim/cv_clim_era5.nc' # '/home/christoph/phd/data/framework_example_ds/aew/' # '/project/meteo/w2w/C3/fischer/belanger/aew_clim/' # -in_files = '/home/ws/he7273/phd_all/data/coll_oper/jja2021/jja2021.nc' # 'C:\\Users\\Christoph\\phd\\data\\enstools-feature\\2008_sum_uv.nc' # '/home/christoph/phd/data/framework_example_ds/aew/cv_aug_08.nc' +cv_clim_file = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/climatologies/cv_climatology_era5/cv_clim_era5.nc' -out_dir = join('/home/ws/he7273/phd_all/data/aew/out/') # '/project/meteo/w2w/C3/fischer/belanger/out/' -plot_dir = join('/home/ws/he7273/phd_all/data/aew/plots/') # '/project/meteo/w2w/C3/fischer/belanger/plots/' - -timedelta_ana = np.timedelta64(7, 'D') - -# aew_kitweather_ecmwf_dir = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/ecmwf-hres/forecasts/' -# 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_' - -# 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_ecmwf = '/home/iconeps/plots/aew/ecmwf-hres/' # '/project/meteo/w2w/C3/fischer/belanger/plots/' plot_dir_ecmwf_prefix = 'aew_ecmwf-hres_' +plot_dir_aifs = '/home/iconeps/plots/aew/ecmwf-aifs/' # '/project/meteo/w2w/C3/fischer/belanger/plots/' +plot_dir_aifs_prefix = 'aew_ecmwf-aifs_' plot_dir_ecmwf_ens = '/home/iconeps/plots/aew/ecmwf-ens/' # '/project/meteo/w2w/C3/fischer/belanger/plots/' plot_dir_ecmwf_ens_prefix = 'aew_ecmwf-ens_' plot_dir_icon = '/home/iconeps/plots/aew/icon-global-det/' # '/project/meteo/w2w/C3/fischer/belanger/plots/' @@ -42,17 +29,18 @@ plot_dir_gfs_ens_prefix = 'aew_gefs_' ens_rain_threshold = 1.0 # display prob of exceeding X mm hr-1 latest_run_dir = '/home/iconeps/plots/aew/' -latest_run_info_file_name = 'latest_aew_run.txt' +latest_run_det_info_file_name = 'latest_aew_det_run.txt' +latest_run_ens_info_file_name = 'latest_aew_ens_run.txt' # out_dir = join(expanduser("~") + '/phd/data/aew/out/') # '/project/meteo/w2w/C3/fischer/belanger/out/' -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_dir = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/climatologies/raw_cv_era5/' # reference where ALL cv data is (only for clim calc.) cv_data_ex = join(cv_data_dir, '*.nc') def get_clim_file(): - fn = aew_clim_dir + fn = cv_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') + "_" @@ -118,4 +106,4 @@ def is_point_in_area(lon, lat): if 0 < lat < 26: return True - return False \ No newline at end of file + return False diff --git a/enstools/feature/identification/african_easterly_waves/identification.py b/enstools/feature/identification/african_easterly_waves/identification.py index 8bd9b22d88839a5bf9f1ec1d1a0bd850abb37986..a0e19a775dac0afba96bb14ede12b60aae39bcee 100644 --- a/enstools/feature/identification/african_easterly_waves/identification.py +++ b/enstools/feature/identification/african_easterly_waves/identification.py @@ -234,6 +234,7 @@ class AEWIdentification(IdentificationStrategy): subplot_kws={'projection': ccrs.PlateCarree()}) paths = c.collections[0].get_paths() + plt.close(fig) wt = data_chunk.wavetroughs diff --git a/enstools/feature/identification/african_easterly_waves/plotting.py b/enstools/feature/identification/african_easterly_waves/plotting.py index 8a49b0667eaad0c56d37d2784a923484fdc0c951..b7191211433dbe9c37a80a51c3e0c4f7b462cc05 100644 --- a/enstools/feature/identification/african_easterly_waves/plotting.py +++ b/enstools/feature/identification/african_easterly_waves/plotting.py @@ -6,6 +6,7 @@ from PIL import Image import matplotlib import matplotlib as mpl import matplotlib.patches as patches +from cartopy.util import add_cyclic_point import cartopy.crs as ccrs import cartopy.feature as cfeature from datetime import datetime @@ -16,7 +17,7 @@ from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER def get_kitweather_rain_cm(ens_mode): rgb_colors = [] - pathtotxtfile = '/home/iconeps/icon_data/additional_data/colorpalettes/' + pathtotxtfile = '/home/iconeps/filter_scripts/aew_script/' if ens_mode: filename_colorpalette = 'colormap_WhiteBeigeGreenBlue_16.txt' else: @@ -31,7 +32,7 @@ def get_kitweather_rain_cm(ens_mode): cmap = cmap.with_extremes(bad='white', under=rgb_colors[0], over=rgb_colors[-1]) if ens_mode: - levels = [0.01,0.02,0.05,0.1,0.12,0.15,0.2,0.25,0.3,0.4,0.5,0.75] + levels = [0.01,0.02,0.05,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8,1.0] else: levels = [0.1,0.2,0.3,0.5,1,2,3,5,10,20,30,50] norm = mpl.colors.BoundaryNorm(levels, cmap.N) @@ -98,14 +99,19 @@ def plot_ts_filtered_waves(wts_part_of_tracks, fig_name, ds=None, tp=None, ens_m ds.plot.streamplot(x='lon', y='lat', u='u', v='v', linewidth=0.6, arrowsize = 0.5, density=6, - color='black') # , transform_first=True not working, or is already implemented. still slow. + color='black') # not working, or is already implemented. still slow. t3 = timer() if tp is not None: - # transform to mm - tp.plot.contourf(levels=levels_rain, extend='max', subplot_kws={'transform_first': True}, - cmap=rain_cm, norm=norm, add_colorbar=False) + tp_wrap_data, tp_wrap_lon = add_cyclic_point(tp.data, coord=tp.lon) + + X, Y = np.meshgrid(tp_wrap_lon, tp.lat) + # use plt API to use wrapped latlonlev data + plt.contourf(X, Y, tp_wrap_data, levels=levels_rain, extend='max', cmap=rain_cm, norm=norm, transform_first=True) + + # tp.plot.contourf(levels=levels_rain, extend='max', subplot_kws={'transform_first': True}, + # cmap=rain_cm, norm=norm, add_colorbar=False) t4 = timer() @@ -142,24 +148,18 @@ def plot_ts_filtered_waves(wts_part_of_tracks, fig_name, ds=None, tp=None, ens_m plt.savefig(fig_name, format='png', backend='agg') - plt.figure().clear() - plt.close() - plt.cla() - plt.clf() + plt.close(fig) crop_top_bottom_whitespace(fig_name) t6 = timer() - """ - print("Init: " + str(t2 - t1)) - print("Streamplot: " + str(t3 - t2)) - print("Rain: " + str(t4 - t3)) - print("Wavetroughs: " + str(t5 - t4)) - print("Finalize: " + str(t6 - t5)) + + #print("Init: " + str(t2 - t1)) + #print("Streamplot: " + str(t3 - t2)) + #print("Rain: " + str(t4 - t3)) + #print("Wavetroughs: " + str(t5 - t4)) + #print("Finalize: " + str(t6 - t5)) - print("Saved to " + fig_name) - exit() - """ print("Saved to " + fig_name) return @@ -461,8 +461,10 @@ def plot_kw(tracks, ds, tp=None, plot_prefix=None, ens_mode=False): os.makedirs(plot_dir, exist_ok=True) if ens_mode and tp is not None: + print("Computing mean precip...") tp_thr = (tp > cfg.ens_rain_threshold).astype(dtype=float) - tp_prob = tp_thr.mean(dim="member") # TODO get_member_dim? + tp_prob = tp_thr.mean(dim="member").compute() # TODO get_member_dim? + print("Done.") tp = tp_prob # call plotting for each date @@ -476,6 +478,7 @@ def plot_kw(tracks, ds, tp=None, plot_prefix=None, ens_mode=False): except (KeyError, AttributeError) as e: print("No rain data for " + str(date)) tp_ss = None + plot_ts_filtered_waves(wavetroughs_in_tracks[date], fig_name, ds=ds_ss, tp=tp_ss, ens_mode=ens_mode) diff --git a/enstools/feature/identification/african_easterly_waves/run_identify.py b/enstools/feature/identification/african_easterly_waves/run_identify.py index 2528313c5244a373977da6f5ecc7aae5db430059..baf343972ef340750f27f9c4227fa68854190554 100644 --- a/enstools/feature/identification/african_easterly_waves/run_identify.py +++ b/enstools/feature/identification/african_easterly_waves/run_identify.py @@ -35,6 +35,11 @@ if len(sys.argv) >= 3 and sys.argv[1] == '-kw': plot_dir = cfg.plot_dir_ecmwf plot_file_prefix = cfg.plot_dir_ecmwf_prefix ens_mode = False + elif sys.argv[2] == 'aifs_fc': + ds, run = get_aifs_forecast(run=run, include_analysis_delta=cfg.timedelta_ana, add_prec_rate=True) + plot_dir = cfg.plot_dir_aifs + plot_file_prefix = cfg.plot_dir_aifs_prefix + ens_mode = False elif sys.argv[2] == 'icon_fc': ds, run = get_icon_forecast(run=run, include_analysis_delta=cfg.timedelta_ana, add_prec_rate=True) plot_dir = cfg.plot_dir_icon @@ -47,11 +52,15 @@ if len(sys.argv) >= 3 and sys.argv[1] == '-kw': ens_mode = False elif sys.argv[2] == 'ecmwf_ens': ds, run = get_ecmwf_ensemble(run=run, add_prec_rate=True) + # remove init time, has no rainfall in some models + ds = [ds_.isel(time=slice(1,None)) for ds_ in ds] plot_dir = cfg.plot_dir_ecmwf_ens plot_file_prefix = cfg.plot_dir_ecmwf_ens_prefix ens_mode = True elif sys.argv[2] == 'gfs_ens': ds, run = get_gfs_ensemble(run=run, add_prec_rate=True) + # remove init time, has no rainfall in some models + ds[0] = ds[0].isel(time=slice(1,None)) # GEFS has no 000h rain, cut first ds0 plot_dir = cfg.plot_dir_gfs_ens plot_file_prefix = cfg.plot_dir_gfs_ens_prefix ens_mode = True @@ -121,13 +130,20 @@ if len(sys.argv) > 1 and sys.argv[1] == '-kw': # All done. Update text file containing time of latest finished run. - if os.path.exists(cfg.plot_dir_ecmwf + run) and os.path.exists(cfg.plot_dir_icon + run) and os.path.exists(cfg.plot_dir_gfs + run) and not ens_mode: - print("Run " + run + " done for ECMWF, ICON, GFS. Update latest run info file.") - yyyymmddhh = run[4:] + yyyymmddhh = run[4:] + if ens_mode and os.path.exists(cfg.plot_dir_ecmwf_ens + run): + print("Ensemble run " + run + " done for ECMWF, GFS. Update latest run info file.") + latest_run_ens_file_path = cfg.latest_run_dir + cfg.latest_run_ens_info_file_name + + with open(latest_run_ens_file_path, 'w+') as info_file: + info_file.write(yyyymmddhh) + - latest_run_file_path = cfg.latest_run_dir + cfg.latest_run_info_file_name + if (not ens_mode) and (os.path.exists(cfg.plot_dir_ecmwf + run) or os.path.exists(cfg.plot_dir_icon + run) or os.path.exists(cfg.plot_dir_gfs + run)): + print("Det run " + run + " done for ECMWF, ICON, GFS. Update latest run info file.") + latest_run_det_file_path = cfg.latest_run_dir + cfg.latest_run_det_info_file_name - with open(latest_run_file_path, 'w+') as info_file: + with open(latest_run_det_file_path, 'w+') as info_file: info_file.write(yyyymmddhh) # All done. scp data over to webserver. diff --git a/enstools/feature/tracking/african_easterly_waves/tracking.py b/enstools/feature/tracking/african_easterly_waves/tracking.py index af1c6439da23384fb20acc91782da2465e89b832..9bd192621d303234366b3f345ae2e5986f1896f7 100644 --- a/enstools/feature/tracking/african_easterly_waves/tracking.py +++ b/enstools/feature/tracking/african_easterly_waves/tracking.py @@ -128,7 +128,8 @@ class AEWTracking(ObjectComparisonTracking): continue # if nodes along track are predominantly horizontal orientated discard them. - lon_lat_ratio = [extent_lons[i] / extent_lats[i] for i in way_node_idxs] + + lon_lat_ratio = [extent_lons[i] / extent_lats[i] if extent_lats[i] > 0.1 else extent_lons[i] for i in way_node_idxs] if statistics.mean(lon_lat_ratio) > 1: continue