From f0a44cfbb6d6106f989eae241a3f8618aa9a150b Mon Sep 17 00:00:00 2001 From: Morten Kretschmer Date: Sun, 21 Apr 2024 13:38:31 +0200 Subject: [PATCH] up-to-date --- .gitignore | 7 +- src/download_script.sh | 4 +- src/main.py | 97 ++++++++++++++----------- src/modelinfolib.py | 89 ++++++++++------------- src/plotlib.py | 161 ++++++++++++++++++++++------------------- src/run.yml | 2 +- src/utilitylib.py | 63 +++++++++------- 7 files changed, 225 insertions(+), 198 deletions(-) diff --git a/.gitignore b/.gitignore index 96cb0ae..16614be 100644 --- a/.gitignore +++ b/.gitignore @@ -1,7 +1,8 @@ # ignore ALL .log files -.ipynb_checkpoints/ -*/.ipynb_checkpoints/ +**__pycache__ + +**.ipynb_checkpoints/ */*.png */*/*.png @@ -17,4 +18,4 @@ */*/*.pyc src/invariant -src/iconnest +src/modeldata diff --git a/src/download_script.sh b/src/download_script.sh index 116fc04..42dcb5c 100755 --- a/src/download_script.sh +++ b/src/download_script.sh @@ -4,7 +4,7 @@ source /etc/profile conda activate HodographMaps -store_path=$(pwd)/iconnest +store_path=$(pwd)/modeldata # create nwp directory and if not there a output images directory mkdir -p ${store_path} @@ -65,4 +65,4 @@ do done # remove nwp files -rm -rf ${store_path} +# rm -rf ${store_path} diff --git a/src/main.py b/src/main.py index 292504a..2e36250 100644 --- a/src/main.py +++ b/src/main.py @@ -11,29 +11,28 @@ import utilitylib as ut import plotlib import modelinfolib as model -model = model.icon_nest # --------------------------------------------------------------------------------------------------------------------- -def run(program_mode, fieldname, rundate, model_run, fp): +def run(model_obj, program_mode, fieldname, fp): config = ut.load_yaml('config.yml') if program_mode == "Test": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./iconnest/") - assert cape_fld.shape == (model.getnlat, model.getnlon), "Shape inconsistency" - plotlib.test_plot(cape_fld, lats, lons, fp, model_run, titel='CAPE') + cape_fld, lats, lons = ut.open_gribfile_single(model_obj, fieldname, fp, path="./modeldata/") + assert cape_fld.shape == (model_obj.getnlat(), model_obj.getnlon()), "Shape inconsistency" + plotlib.test_plot(model_obj, cape_fld, lats, lons, fp) elif program_mode == "Sounding": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./iconnest/") + cape_fld, lats, lons = ut.open_gribfile_single(model_obj, fieldname, fp, path="./modeldata/") steps = config["steps"] - nlvl = int(model.getnlev()/steps) - t_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + nlvl = int(model_obj.getnlev()/steps) + t_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) t_fld.fill(np.nan) - q_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + q_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) q_fld.fill(np.nan) - p_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + p_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) p_fld.fill(np.nan) if config["levels"][0] > config["levels"][1]: @@ -44,24 +43,24 @@ def run(program_mode, fieldname, rundate, model_run, fp): lvl_idx = 0 for level in range(config["levels"][0], config["levels"][1], steps): - t_fld[lvl_idx, :, :] = ut.open_gribfile_multi("T", level, rundate, model_run, fp, path="./iconnest/") - q_fld[lvl_idx, :, :] = ut.open_gribfile_multi("QV", level, rundate, model_run, fp, path="./iconnest/") - p_fld[lvl_idx, :, :] = ut.open_gribfile_multi("P", level, rundate, model_run, fp, path="./iconnest/") + t_fld[lvl_idx, :, :] = ut.open_gribfile_multi(model_obj, "T", level, fp, path="./modeldata/") + q_fld[lvl_idx, :, :] = ut.open_gribfile_multi(model_obj, "QV", level, fp, path="./modeldata/") + p_fld[lvl_idx, :, :] = ut.open_gribfile_multi(model_obj, "P", level, fp, path="./modeldata/") lvl_idx += 1 if lvl_idx >= nlvl: break print(np.nanmean(t_fld, axis=(1, 2))-273.15) - plotlib.sounding_plot(cape_fld, t_fld, q_fld, p_fld, lats, lons, fp, model_run, titel='CAPE') + plotlib.sounding_plot(model_obj, cape_fld, t_fld, q_fld, p_fld, lats, lons, fp) elif program_mode == "Basic": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./iconnest/") + cape_fld, lats, lons = ut.open_gribfile_single(model_obj, fieldname, fp, path="./modeldata/") steps = config["steps"] - nlvl = int(model.getnlev()/steps) - u_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + nlvl = int(model_obj.getnlev()/steps) + u_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) u_fld.fill(np.nan) - v_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + v_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) v_fld.fill(np.nan) if config["levels"][0] > config["levels"][1]: @@ -72,30 +71,30 @@ def run(program_mode, fieldname, rundate, model_run, fp): lvl_idx = 0 for level in range(config["levels"][0], config["levels"][1], steps): - u_fld[lvl_idx, :, :] = ut.open_gribfile_multi("U", level, rundate, model_run, fp, path="./iconnest/") - v_fld[lvl_idx, :, :] = ut.open_gribfile_multi("V", level, rundate, model_run, fp, path="./iconnest/") + u_fld[lvl_idx, :, :] = ut.open_gribfile_multi(model_obj, "U", level, fp, path="./modeldata/") + v_fld[lvl_idx, :, :] = ut.open_gribfile_multi(model_obj, "V", level, fp, path="./modeldata/") lvl_idx += 1 if lvl_idx >= nlvl: break print(np.nanmean(u_fld, axis=(1, 2))) - plotlib.basic_plot(cape_fld, u_fld, v_fld, lats, lons, fp, model_run, - titel='CAPE', threshold=config["threshold"]) - plotlib.basic_plot_custarea(cape_fld, u_fld, v_fld, lats, lons, fp, model_run, - titel='CAPE', threshold=config["threshold"]) + plotlib.basic_plot(model_obj, cape_fld, u_fld, v_fld, lats, lons, fp, + threshold=config["threshold"]) + plotlib.basic_plot_custarea(model_obj, cape_fld, u_fld, v_fld, lats, lons, fp, + threshold=config["threshold"]) elif program_mode == "Nixon": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./iconnest/") + cape_fld, lats, lons = ut.open_gribfile_single(model_obj, fieldname, fp, path="./modeldata/") steps = config["steps"] - nlvl = int(model.getnlev()/steps) - u_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + nlvl = int(model_obj.getnlev()/steps) + u_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) u_fld.fill(np.nan) - v_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + v_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) v_fld.fill(np.nan) - h_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + h_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) h_fld.fill(np.nan) - p_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + p_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) p_fld.fill(np.nan) if config["levels"][0] > config["levels"][1]: @@ -106,10 +105,10 @@ def run(program_mode, fieldname, rundate, model_run, fp): lvl_idx = 0 for level in range(config["levels"][0], config["levels"][1], steps): - u_fld[lvl_idx, :, :] = ut.open_gribfile_multi("U", level, rundate, model_run, fp, path="./iconnest/") - v_fld[lvl_idx, :, :] = ut.open_gribfile_multi("V", level, rundate, model_run, fp, path="./iconnest/") - h_fld[lvl_idx, :, :] = ut.open_gribfile_multi("H", level, rundate, model_run, fp, path="./iconnest/") - p_fld[lvl_idx, :, :] = ut.open_gribfile_multi("P", level, rundate, model_run, fp, path="./iconnest/") + u_fld[lvl_idx, :, :] = ut.open_gribfile_multi(model_obj, "U", level, fp, path="./modeldata/") + v_fld[lvl_idx, :, :] = ut.open_gribfile_multi(model_obj, "V", level, fp, path="./modeldata/") + h_fld[lvl_idx, :, :] = ut.open_gribfile_multi(model_obj, "H", level, fp, path="./modeldata/") + p_fld[lvl_idx, :, :] = ut.open_gribfile_multi(model_obj, "P", level, fp, path="./modeldata/") lvl_idx += 1 if lvl_idx >= nlvl: @@ -120,7 +119,7 @@ def run(program_mode, fieldname, rundate, model_run, fp): du = np.subtract(u_fld[30, :, :], u_fld[0, :, :]) dv = np.subtract(v_fld[30, :, :], v_fld[0, :, :]) dls_fld = np.sqrt(np.add(np.square(du), np.square(dv))) - plotlib.nixon_proj(cape_fld, dls_fld, u_fld, v_fld, p_fld, h_fld, lats, lons, fp, model_run, imfmt="png") + plotlib.nixon_proj(model_obj, cape_fld, dls_fld, u_fld, v_fld, p_fld, h_fld, lats, lons, fp, imfmt="png") else: print("Wrong command line argument") exit(-1) @@ -151,6 +150,9 @@ def main(): parser.add_argument('-r', '--run', type=str, help='Run (0,6,12,18,.. .etc)') + parser.add_argument('-m', '--model', type=str, + help='--model: ICON, IFS, or GFS') + # Parse the command line arguments args = parser.parse_args() @@ -161,10 +163,22 @@ def main(): print("Unknown field") exit(-1) + if args.model is None: + model_obj = model.MODELIFNO("ICON EU", 1377, 657, 0.0625, 54, "pres") + elif "ICON" in args.model: + model_obj = model.icon_nest + # elif args.model == "IFS": + # model_obj = model.ifs + # elif args.model == "GFS": + # model_obj = model.gfs + else: + print("Unkown model! Exit.") + exit(0) + if args.date is None: - rundate = datetime.strptime(run_config["default_date"], "%Y-%m-%d") + model_obj.setrundate(datetime.strptime(run_config["default_date"], "%Y-%m-%d")) else: - rundate = datetime.strptime(args.date, "%Y-%m-%d") + model_obj.setrundate(datetime.strptime(args.date, "%Y-%m-%d")) if args.fp is None: fp = run_config["fp"] @@ -172,9 +186,9 @@ def main(): fp = args.fp if args.run is None: - model_run = run_config["run"] + model_obj.setrun(run_config["run"]) else: - model_run = args.run + model_obj.setrun(args.run) if args.mode != "Test" and args.mode != "Sounding" and args.mode != "Basic" and args.mode != "Nixon": print("Unknown Mode. Exit program.") @@ -185,9 +199,8 @@ def main(): # replace space with underscores fieldname = args.field.replace(" ", "_") - print(f"\nDate: {rundate}\n Arguments: {args} \nConfig-File: {config}\n\n") - - run(program_mode, fieldname, rundate, model_run, fp) + print(f"\nArguments: {args}\n{model_obj}") + run(model_obj, program_mode, fieldname, fp) # --------------------------------------------------------------------------------------------------------------------- diff --git a/src/modelinfolib.py b/src/modelinfolib.py index 67e3402..ac1bafc 100644 --- a/src/modelinfolib.py +++ b/src/modelinfolib.py @@ -1,61 +1,22 @@ #!/usr/bin/python3 -""" -# ICON Nest -points = 904689 -nlon = 1377 -nlat = 657 -nlev = 51 -lonmin = -23.5 -lonmax = 45.0 -latmin = 29.5 -latmax = 70.5 -d_grad = 0.0625 - -# IFS -points = 405900 -nlon = 900 -nlat = 451 -nlev = -lonmin = 0 -lonmax = 359 -latmin = -90 -latmax = 90 -d_grad = 0.4 - -# GFS -points = 1038240 -nlon = 1440 -nlat = 721 -nlev = -lonmin = 0 -lonmax = 359 -latmin = -90 -latmax = 90 -d_grad = 0.25 - -""" +from datetime import datetime, date class MODELIFNO: - def __init__(self, nlon, nlat, d_grad, nlev, levtyp): + def __init__(self, modelname, nlon, nlat, d_grad, nlev, levtyp): + self.modelname = modelname self.points = nlon*nlat self.nlon = nlon self.nlat = nlat self.nlev = nlev self.d_grad = d_grad self.levtyp = levtyp - - """ - def __init__(self, nlon, nlat, nlev, lonmin, lonmax, latmin, latmax, d_grad): - self.points = nlon*nlat - self.nlon = nlon - self.nlat = nlat - self.nlev = nlev - self.d_grad = d_grad - print(lonmin, lonmax, latmin, latmax) - """ - + self.run = -99 + self.rundate = date.today() + self.hodo_interval_lat = range(272, 415, 12) + self.hodo_interval_lon = range(420, 670, 15) + def __str__(self): return ( f"Model Information:\nPoints: {self.points}\n" @@ -64,6 +25,33 @@ def __str__(self): f"Number of Levels: {self.nlev}\tLeveltyps: {self.levtyp}\n" ) + def setrun(self, run): + if isinstance(run, int): + self.run = run + else: + self.run = int(run) + + def setrundate(self, rundate): + if isinstance(rundate, datetime): + self.rundate = rundate + else: + self.rundate = datetime.strptime(rundate, "%Y-%m-%d") + + def getrun(self): + return self.run + + def getrundate(self): + return self.rundate + + def getrundate_as_str(self, fmt="%Y%m%d"): + return self.rundate.strftime(fmt) + + def getname(self): + return self.modelname + + def getParamter(self): + return self.parlist + def getpoints(self): return self.points @@ -74,7 +62,7 @@ def getnlat(self): return self.nlat def getnlev(self): - return self.nlev + return int(self.nlev) def getlevtyp(self): return self.d_grad @@ -82,6 +70,9 @@ def getlevtyp(self): def getd_grad(self): return self.d_grad + def create_plottitle(self): + return f"Hodographmap of {self.modelname}" + # Example usage: -icon_nest = MODELIFNO(1377, 657, 0.0625, 54, "model") # lowest 74 and we download till 20 +icon_nest = MODELIFNO("ICON Nest", 1377, 657, 0.0625, 54, "model") # lowest 74 and we download till 20 diff --git a/src/plotlib.py b/src/plotlib.py index 60351a3..b076092 100644 --- a/src/plotlib.py +++ b/src/plotlib.py @@ -33,7 +33,7 @@ titlesize = config["titlesize"] -def eu_merc(hour, start, projection=crs.Mercator(), factor=3): +def eu_merc(hour, start, datetime_obj, model_name, projection=crs.Mercator(), factor=3): fig, ax = plt.subplots(figsize=(3*factor, 3.5091*factor), subplot_kw=dict(projection=projection)) ax.set_extent([-10.5, 28.0, 30.5, 67.5]) # ax.stock_img() @@ -43,60 +43,59 @@ def eu_merc(hour, start, projection=crs.Mercator(), factor=3): gl.ylabels_right = False gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER - string1, string2 = ut.datum(hour, start) - plt.annotate(string1, xy=(0.8, 1), xycoords='axes fraction', fontsize=fontsize) - plt.annotate(string2, xy=(0, 1), xycoords='axes fraction', fontsize=fontsize) + string1, string2 = ut.datum(hour, start, datetime_obj) + plt.annotate(model_name, xy=(0.7, -0.04), xycoords='axes fraction', fontsize=10) + plt.annotate(string1, xy=(0.82, 1.01), xycoords='axes fraction', fontsize=fontsize) + plt.annotate(string2, xy=(0, 1.01), xycoords='axes fraction', fontsize=fontsize) return fig, ax -def eu_states(hour, start, projection=crs.EuroPP()): +def eu_states(hour, start, datetime_obj, model_name, projection=crs.EuroPP()): fig, ax = plt.subplots(figsize=(9, 7), subplot_kw=dict(projection=projection)) plt.subplots_adjust(left=0.05, right=0.97, bottom=0.05, top=0.95) ax.set_extent([-9.5, 33.0, 38.5, 58.5]) ax.coastlines('50m', linewidth=1.2) ax.add_feature(states_provinces, edgecolor='black') - string1, string2 = ut.datum(hour, start) - plt.annotate("ICON Nest (DWD)", xy=(0.7, -0.04), xycoords='axes fraction', fontsize=10) - plt.annotate(string1, xy=(0.8, 1), xycoords='axes fraction', fontsize=fontsize) - plt.annotate(string2, xy=(0, 1), xycoords='axes fraction', fontsize=fontsize) + string1, string2 = ut.datum(hour, start, datetime_obj) + plt.annotate(model_name, xy=(0.7, -0.04), xycoords='axes fraction', fontsize=10) + plt.annotate(string1, xy=(0.82, 1.01), xycoords='axes fraction', fontsize=fontsize) + plt.annotate(string2, xy=(0, 1.01), xycoords='axes fraction', fontsize=fontsize) return fig, ax -def ce_states(hour, start, projection=crs.EuroPP(), lon1=1.56, lon2=18.5, lat1=45.1, lat2=56.6): +def ce_states(hour, start, datetime_obj, projection=crs.EuroPP(), lon1=1.56, lon2=18.5, lat1=45.1, lat2=56.6): fig, ax = plt.subplots(figsize=(11, 9), subplot_kw=dict(projection=projection)) plt.subplots_adjust(left=0.05, right=0.99, bottom=0.05, top=0.95) ax.set_extent([lon1, lon2, lat1, lat2]) ax.coastlines('10m', linewidth=1.2) ax.add_feature(states_provinces, edgecolor='black') - string1, string2 = ut.datum(hour, start) - plt.annotate("ICON Nest (DWD)", xy=(0.02, -0.04), xycoords='axes fraction', fontsize=10) - plt.annotate(string1, xy=(0.835, 1.02), xycoords='axes fraction', fontsize=fontsize) - plt.annotate(string2, xy=(0, 1.02), xycoords='axes fraction', fontsize=fontsize) + string1, string2 = ut.datum(hour, start, datetime_obj) + plt.annotate(string1, xy=(0.82, 1.01), xycoords='axes fraction', fontsize=fontsize) + plt.annotate(string2, xy=(0, 1.01), xycoords='axes fraction', fontsize=fontsize) return fig, ax -def customize_area(hour, start, projection=crs.EuroPP(), lon1=10.7, lon2=18, lat1=49.8, lat2=54.8): +def customize_area(hour, start, datetime_obj, projection=crs.EuroPP(), lon1=10.7, lon2=18, lat1=49.8, lat2=54.8): fig, ax = plt.subplots(figsize=(11, 9), subplot_kw=dict(projection=projection)) plt.subplots_adjust(left=0.05, right=0.99, bottom=0.1, top=0.95) ax.set_extent([lon1, lon2, lat1, lat2]) ax.coastlines('10m', linewidth=1.2) ax.add_feature(states_provinces, edgecolor='black') - string1, string2 = ut.datum(hour, start) - plt.annotate("ICON Nest (DWD)", xy=(0.02, -0.04), xycoords='axes fraction', fontsize=10) + string1, string2 = ut.datum(hour, start, datetime_obj) plt.annotate(string1, xy=(0.82, 1.01), xycoords='axes fraction', fontsize=fontsize) plt.annotate(string2, xy=(0, 1.01), xycoords='axes fraction', fontsize=fontsize) return fig, ax -def alps(hour, start, projection=crs.EuroPP(), lon1=5.8, lon2=17.8, lat1=45.23, lat2=49.5): +def alps(hour, start, datetime_obj, projection=crs.EuroPP(), lon1=5.8, lon2=17.8, lat1=45.23, lat2=49.5): fig, ax = plt.subplots(figsize=(11, 9), subplot_kw=dict(projection=projection)) ax.set_extent([lon1, lon2, lat1, lat2]) ax.coastlines('10m', linewidth=1.2) ax.add_feature(states_provinces, edgecolor='black') - string1, string2 = ut.datum(hour, start) + string1, string2 = ut.datum(hour, start, datetime_obj) plt.annotate("ICON Nest (DWD)", xy=(0.7, -0.04), xycoords='axes fraction', fontsize=10) - plt.annotate(string1, xy=(0.8, 1), xycoords='axes fraction', fontsize=fontsize) - plt.annotate(string2, xy=(0, 1), xycoords='axes fraction', fontsize=fontsize) + plt.annotate(string1, xy=(0.82, 1.01), xycoords='axes fraction', fontsize=fontsize) + plt.annotate(string2, xy=(0, 1.01), xycoords='axes fraction', fontsize=fontsize) return fig, ax @@ -123,7 +122,7 @@ def two_plots(projection=crs.EuroPP(), lon1=3.56, lon2=16.5, lat1=46.2, lat2=55. # --------------------------------------------------------------------------------------------------------------------- -def test_plot(cape_fld, lats, lons, hour, run, titel='CAPE'): +def test_plot(model_obj, cape_fld, lats, lons, hour): """ Parameters: ------------ @@ -137,9 +136,11 @@ def test_plot(cape_fld, lats, lons, hour, run, titel='CAPE'): -------- None """ + model_name = model_obj.getname() + start = model_obj.getrun() - fig, ax = ce_states(hour, run, projection=crs.PlateCarree()) - plt.title(titel, fontsize=titlesize) + fig, ax = ce_states(hour, start, model_obj.getrundate(), projection=crs.PlateCarree()) + plt.title(model_obj.create_plottitle(), fontsize=titlesize) wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), cmap=cmap, extend='max', alpha=0.4, antialiased=True) @@ -147,7 +148,7 @@ def test_plot(cape_fld, lats, lons, hour, run, titel='CAPE'): cax = fig.add_axes([0.27, 0.05, 0.35, 0.05]) fig.colorbar(wx, cax=cax, orientation='horizontal') - name = f"test_ce_{hour}.png" + name = f"test_{model_name}_{hour}.png" plt.savefig(name) plt.close() @@ -207,7 +208,7 @@ def soundingpoint(point, temp, q, pres, ax, width=0.1, proj='skewT', smooth=Fals # --------------------------------------------------------------------------------------------------------------------- -def sounding_plot(cape_fld, temp, q_fld, p_fld, lats, lons, hour, start, titel='CAPE', imfmt="png"): +def sounding_plot(model_obj, cape_fld, temp, q_fld, p_fld, lats, lons, hour, imfmt="png"): """ Parameters: ------------ @@ -223,8 +224,11 @@ def sounding_plot(cape_fld, temp, q_fld, p_fld, lats, lons, hour, start, titel=' None """ - fig, ax = ce_states(hour, start, projection=crs.PlateCarree()) - plt.title(titel, fontsize=titlesize) + model_name = model_obj.getname() + start = model_obj.getrun() + + fig, ax = ce_states(hour, start, model_obj.getrundate(), projection=crs.PlateCarree()) + plt.title(model_obj.create_plottitle(), fontsize=titlesize) wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), cmap=cmap2, extend='max', alpha=0.4, antialiased=True) @@ -235,7 +239,7 @@ def sounding_plot(cape_fld, temp, q_fld, p_fld, lats, lons, hour, start, titel=' np.mean(temp[:, i-1:i+1, j-1:j+1], axis=(1, 2)), np.mean(q_fld[:, i-1:i+1, j-1:j+1], axis=(1, 2)), np.mean(p_fld[:, i-1:i+1, j-1:j+1], axis=(1, 2)), - ax, width=0.05) # , proj=crs.PlateCarree() + ax, width=0.05) cax = fig.add_axes([0.27, 0.05, 0.35, 0.05]) fig.colorbar(wx, cax=cax, orientation='horizontal') @@ -243,7 +247,7 @@ def sounding_plot(cape_fld, temp, q_fld, p_fld, lats, lons, hour, start, titel=' ax.annotate('red dot line: freezing level', xy=(0.75, -0.1), xycoords='axes fraction', fontsize=14) ax.annotate('blue lines: temperature', xy=(0.75, -0.04), xycoords='axes fraction', fontsize=14) ax.annotate('green lines: dewpoints', xy=(0.75, -0.07), xycoords='axes fraction', fontsize=14) - name = f"./images/soundingmap_ce_{hour}.{imfmt}" + name = f"./images/soundingmap_{model_name}_{hour}.{imfmt}" plt.savefig(name) plt.close() @@ -286,7 +290,6 @@ def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): ax2.plot(np.linspace(0, 2*np.pi, 100), np.zeros(100)+10, '-k', alpha=.3, lw=0.8) # ax2.plot(np.linspace(0, 2*np.pi, 100), np.zeros(100)+30, '-k', alpha=.3, lw=0.8) - # plot data wdir, spd = met.uv2spddir(u, v) # smoothing @@ -294,10 +297,14 @@ def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): wdir[1:] = (wdir[1:] + wdir[:-1])/2 spd[1:] = (spd[1:] + spd[:-1])/2 + # plot data # draw part of second cricle if np.max(spd[:-20]) > 28: - ax2.plot(np.linspace(np.mean(wdir[np.where(spd[:-20] > 25)])-np.pi/8, - np.mean(wdir[np.where(spd[:-20] > 25)])+np.pi/8, 100), + u_mean = np.mean(u[np.where(spd > 25)]) + v_mean = np.mean(v[np.where(spd > 25)]) + wdir_mean = met.uv2spddir(u_mean, v_mean)[0] + ax2.plot(np.linspace(wdir_mean-np.pi/8, + wdir_mean+np.pi/8, 100), np.zeros(100)+30, '-k', alpha=.3, lw=0.8) ax2.plot(wdir[:10:1], spd[:10:1], 'r-', lw=1.5) @@ -308,7 +315,7 @@ def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): # --------------------------------------------------------------------------------------------------------------------- -def basic_plot(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold=10., imfmt="png"): +def basic_plot(model_obj, cape_fld, u, v, lats, lons, hour, threshold=10., imfmt="png"): """ Parameters: ------------ @@ -323,44 +330,50 @@ def basic_plot(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold= -------- None """ + model_name = model_obj.getname() + start = model_obj.getrun() - fig, ax = ce_states(hour, start, projection=crs.PlateCarree()) - plt.title(titel, fontsize=titlesize) + fig, ax = ce_states(hour, start, model_obj.getrundate(), projection=crs.PlateCarree()) + plt.title(model_obj.create_plottitle(), fontsize=titlesize) wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), cmap=cmap, extend='max', alpha=0.4, antialiased=True) - for i in range(275, 415, 10): - for j in range(420, 670, 15): + for i in model_obj.hodo_interval_lat: + for j in model_obj.hodo_interval_lon: if np.mean(cape_fld[i-1:i+1, j-1:j+1]) > threshold: hodopoint((lons[i, j], lats[i, j]), np.mean(u[:, i-1:i+1, j-1:j+1], axis=(1, 2)), - np.mean(v[:, i-1:i+1, j-1:j+1], axis=(1, 2)), ax, width=0.1) # , proj=crs.PlateCarree() + np.mean(v[:, i-1:i+1, j-1:j+1], axis=(1, 2)), ax, width=0.1) - cax = fig.add_axes([0.27, 0.05, 0.35, 0.05]) + cax = fig.add_axes([0.44, 0.03, 0.35, 0.05]) fig.colorbar(wx, cax=cax, orientation='horizontal') - if "CAPE" in titel: - ax.annotate(r'$J/kg$', xy=(0.65, -0.04), xycoords='axes fraction', fontsize=14) - else: - ax.annotate(r'$m^2/s^2$', xy=(0.65, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('red: 1-10 model level', xy=(0.75, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('green: 10-20 model level', xy=(0.75, -0.07), xycoords='axes fraction', fontsize=14) - ax.annotate('blue: 20-50 model level', xy=(0.75, -0.1), xycoords='axes fraction', fontsize=14) - ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.07), xycoords='axes fraction', fontsize=10) + ax.annotate("CAPE ML(contour plot)", xy=(0.8, -0.07), xycoords='axes fraction', fontsize=14) + ax.annotate(r'in $J/kg$', xy=(0.8, -0.1), xycoords='axes fraction', fontsize=14) + + ax.annotate('red: 1-10 model level', xy=(0.02, -0.03), xycoords='axes fraction', fontsize=14) + ax.annotate('green: 10-20 model level', xy=(0.02, -0.06), xycoords='axes fraction', fontsize=14) + ax.annotate('blue: 20-50 model level', xy=(0.02, -0.09), xycoords='axes fraction', fontsize=14) + ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.12), xycoords='axes fraction', fontsize=10) - name = f"./images/hodographmap_ce_{hour}.{imfmt}" + name = f"./images/hodographmap_{model_name}_{hour}.{imfmt}" plt.savefig(name) plt.close() -def basic_plot_custarea(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold=10., imfmt="png"): +def basic_plot_custarea(model_obj, cape_fld, u, v, lats, lons, hour, threshold=10., imfmt="png"): lon1 = config["customize"]["lon1"] lon2 = config["customize"]["lon2"] lat1 = config["customize"]["lat1"] lat2 = config["customize"]["lat2"] - fig, ax = customize_area(hour, start, projection=crs.PlateCarree(), lon1=lon1, lon2=lon2, lat1=lat1, lat2=lat2) - plt.title(titel, fontsize=titlesize) + + model_name = model_obj.getname() + start = model_obj.getrun() + + fig, ax = ce_states(hour, start, model_obj.getrundate(), + projection=crs.PlateCarree(), lon1=lon1, lon2=lon2, lat1=lat1, lat2=lat2) + plt.title(model_obj.create_plottitle(), fontsize=titlesize) wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), cmap=cmap, extend='max', alpha=0.4, antialiased=True) @@ -372,19 +385,18 @@ def basic_plot_custarea(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', t np.mean(u[:, i-1:i+1, j-1:j+1], axis=(1, 2)), np.mean(v[:, i-1:i+1, j-1:j+1], axis=(1, 2)), ax, width=0.1) # , proj=crs.PlateCarree() - cax = fig.add_axes([0.27, 0.05, 0.35, 0.05]) + cax = fig.add_axes([0.44, 0.03, 0.35, 0.05]) fig.colorbar(wx, cax=cax, orientation='horizontal') - if "CAPE" in titel: - ax.annotate(r'$J/kg$', xy=(0.65, -0.04), xycoords='axes fraction', fontsize=14) - else: - ax.annotate(r'$m^2/s^2$', xy=(0.65, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('red: 1-10 model level', xy=(0.75, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('green: 10-20 model level', xy=(0.75, -0.07), xycoords='axes fraction', fontsize=14) - ax.annotate('blue: 20-50 model level', xy=(0.75, -0.1), xycoords='axes fraction', fontsize=14) - ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.07), xycoords='axes fraction', fontsize=10) + ax.annotate("CAPE ML(contour plot)", xy=(0.8, -0.07), xycoords='axes fraction', fontsize=14) + ax.annotate(r'in $J/kg$', xy=(0.8, -0.1), xycoords='axes fraction', fontsize=14) + + ax.annotate('red: 1-10 model level', xy=(0.02, -0.03), xycoords='axes fraction', fontsize=14) + ax.annotate('green: 10-20 model level', xy=(0.02, -0.06), xycoords='axes fraction', fontsize=14) + ax.annotate('blue: 20-50 model level', xy=(0.02, -0.09), xycoords='axes fraction', fontsize=14) + ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.12), xycoords='axes fraction', fontsize=10) - name = f"./images/hodographmap_area_{hour}.{imfmt}" + name = f"./images/hodographmap_area_{model_name}_{hour}.{imfmt}" plt.savefig(name) plt.close() @@ -449,7 +461,7 @@ def nixon_hodograph(point, u, v, p, height, ax, width=0.1, clim=40, proj='polar' ax2.arrow(theta, 0, 0, mag, head_width=0.1, head_length=0.1) -def nixon_proj(cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt="png"): +def nixon_proj(model_obj, cape_fld, dls_fld, u, v, p, high, lats, lons, hour, imfmt="png"): """ Nixon projection cape_fld : 2D cape field @@ -462,10 +474,11 @@ def nixon_proj(cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt= only hodographs with more the 10 m/s dls """ - titel = 'CAPE with Hodographs' + model_name = model_obj.getname() + start = model_obj.getrun() - fig, ax = ce_states(hour, start, projection=crs.PlateCarree()) - plt.title(titel, fontsize=titlesize) + fig, ax = ce_states(hour, start, model_obj.getrundate(), projection=crs.PlateCarree()) + plt.title(model_obj.create_plottitle(), fontsize=titlesize) clevs = np.array([20, 250, 500, 750, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500]) # clevs = np.array([20, 50, 100, 250, 500, 750, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500]) @@ -489,15 +502,17 @@ def nixon_proj(cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt= np.mean(p[::-1, i-1:i+1, j-1:j+1], axis=(1, 2)), np.mean(high[::-1, i-1:i+1, j-1:j+1], axis=(1, 2)), ax, width=0.1, proj=crs.PlateCarree()) - cax = fig.add_axes([0.27, 0.05, 0.35, 0.05]) + cax = fig.add_axes([0.44, 0.03, 0.35, 0.05]) fig.colorbar(wx, cax=cax, orientation='horizontal') - ax.annotate(r'$J/kg$', xy=(0.65, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('red: 1-10 model level', xy=(0.75, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('green: 10-20 model level', xy=(0.75, -0.07), xycoords='axes fraction', fontsize=14) - ax.annotate('blue: 20-40 model level', xy=(0.75, -0.1), xycoords='axes fraction', fontsize=14) - ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.07), xycoords='axes fraction', fontsize=10) + ax.annotate("CAPE ML(contour plot)", xy=(0.8, -0.07), xycoords='axes fraction', fontsize=14) + ax.annotate(r'in $J/kg$', xy=(0.8, -0.1), xycoords='axes fraction', fontsize=14) + + ax.annotate('red: 1-10 model level', xy=(0.02, -0.03), xycoords='axes fraction', fontsize=14) + ax.annotate('green: 10-20 model level', xy=(0.02, -0.06), xycoords='axes fraction', fontsize=14) + ax.annotate('blue: 20-50 model level', xy=(0.02, -0.09), xycoords='axes fraction', fontsize=14) + ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.12), xycoords='axes fraction', fontsize=10) - name = f"./images/nixon_ce_{hour}.{imfmt}" + name = f"./images/nixon_{model_name}_{hour}.{imfmt}" plt.savefig(name) plt.close() diff --git a/src/run.yml b/src/run.yml index 019cabf..0ef27e1 100644 --- a/src/run.yml +++ b/src/run.yml @@ -1,3 +1,3 @@ run: 0 fp: 15 -default_date: "2024-04-09" +default_date: "2024-04-21" diff --git a/src/utilitylib.py b/src/utilitylib.py index ee4dae0..a27a67e 100644 --- a/src/utilitylib.py +++ b/src/utilitylib.py @@ -1,6 +1,6 @@ #!/usr/bin/python3 -import datetime +from datetime import datetime, timedelta import requests import yaml import pygrib @@ -9,22 +9,20 @@ # --------------------------------------------------------------------------------------------------------------------- -def datum(h, start): - if not isinstance(h, int): - h = int(h) +def datum(leadtime, start, datetime_obj): + if not isinstance(leadtime, int): + leadtime = int(leadtime) if not isinstance(start, int): start = int(start) - today = datetime.date.today() - today = today.timetuple() - x = datetime.datetime(today.tm_year, today.tm_mon, today.tm_mday, start) - # string = "Forecasttime " + x.strftime("%d.%m. %H:%M") + " UTC " - string1 = "Run: " + x.strftime("%d.%m. %H UTC") - x = datetime.datetime(today.tm_year, today.tm_mon, today.tm_mday, start) + datetime.timedelta(hours=h) - string2 = x.strftime("%A, %d.%m. %H UTC") + today = datetime_obj.timetuple() + x = datetime(today.tm_year, today.tm_mon, today.tm_mday, start) + modelrun_string = "Run: " + x.strftime("%d.%m. %H UTC") + x = datetime(today.tm_year, today.tm_mon, today.tm_mday, start) + timedelta(hours=leadtime) + valid_string = x.strftime("%A, %d.%m. %H UTC") - return string1, string2 + return modelrun_string, valid_string # --------------------------------------------------------------------------------------------------------------------- @@ -58,8 +56,13 @@ def download_nwp(fieldname, datum="20240227", run="00", fp=0, store_path="./"): # --------------------------------------------------------------------------------------------------------------------- -def open_gribfile_single(fieldname, datetime_obj, run, fp, path="./iconnest/"): - date_string = datetime_obj.strftime("%Y%m%d") +def open_gribfile_single(model_obj, fieldname, fp, path="./modeldata/", debug=False): + date_string = model_obj.getrundate_as_str("%Y%m%d") + run = model_obj.getrun() + + shape = (model_obj.getnlat(), model_obj.getnlon()) + cape_fld = np.full(shape, np.nan) + nwp_singlelevel = "icon-eu_europe_regular-lat-lon_single-level" # Open the GRIB file @@ -68,21 +71,25 @@ def open_gribfile_single(fieldname, datetime_obj, run, fp, path="./iconnest/"): # Read specific data from the file first_message = grbs[1] - data = first_message.values + cape_fld = first_message.values # Convert missing values to NaNs - data[data == first_message.missingValue] = np.nan + cape_fld[cape_fld == first_message.missingValue] = np.nan lats, lons = first_message.latlons() grbs.close() - # print("Data shape:", data.shape) - print("Maximum:", np.nanmax(data)) - return data, lats, lons + if debug is True: + # print("Data shape:", data.shape) + print("Maximum:", np.nanmax(cape_fld)) + return cape_fld, lats, lons -def open_gribfile_multi(fieldname, lvl, datetime_obj, run, fp, path="./iconnest/"): - date_string = datetime_obj.strftime("%Y%m%d") +def open_gribfile_multi(model_obj, fieldname, lvl, fp, path="./modeldata/", debug=False): + date_string = model_obj.getrundate_as_str("%Y%m%d") + run = model_obj.getrun() nwp_modellevel = "icon-eu_europe_regular-lat-lon_model-level" + shape = (model_obj.getnlat(), model_obj.getnlon()) + data_fld = np.full(shape, np.nan) # Open the GRIB file grbs = pygrib.open(f"{path}{nwp_modellevel}_{date_string}{run:02d}_{fp:03d}_{lvl}_{fieldname.upper()}.grib2") @@ -90,20 +97,20 @@ def open_gribfile_multi(fieldname, lvl, datetime_obj, run, fp, path="./iconnest/ # Read specific data from the file first_message = grbs[1] - data = first_message.values + data_fld = first_message.values # Convert missing values to NaNs - data[data == first_message.missingValue] = np.nan + data_fld[data_fld == first_message.missingValue] = np.nan grbs.close() - - # print("Data shape:", data.shape) - print("Maximum:", np.nanmax(data)) - return data + if debug is True: + # print("Data shape:", data.shape) + print("Maximum:", np.nanmax(data_fld)) + return data_fld # --------------------------------------------------------------------------------------------------------------------- -def open_netcdf(fieldname, path="./iconnest/"): +def open_netcdf(fieldname, path="./modeldata/"): from netCDF4 import Dataset data = Dataset(f"{path}{fieldname}", 'r')