diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..389fea3 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2024 Morten Kretschmer + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md index a32b4a4..acb7d39 100644 --- a/README.md +++ b/README.md @@ -42,3 +42,7 @@ python3 main.py Sounding ## Datasource - ICON Nest: https://opendata.dwd.de/weather/nwp/icon-eu/ - IFS: https://www.ecmwf.int/en/forecasts/datasets/open-data + +## License + +This project is licensed under the terms of the MIT license. See the [LICENSE](LICENSE) file for details. diff --git a/images/hodographmap_ce_9.png b/images/hodographmap_ce_9.png index 390a52f..2572ad0 100644 Binary files a/images/hodographmap_ce_9.png and b/images/hodographmap_ce_9.png differ diff --git a/src/creategridinfo.py b/src/creategridinfo.py index 3aee04b..c728a5f 100755 --- a/src/creategridinfo.py +++ b/src/creategridinfo.py @@ -3,23 +3,29 @@ import numpy as np -path_icon='./' +path_icon = './' if (os.getcwd() != path_icon): os.chdir(path_icon) + def creategridinfo(d_grad=0.0625, filename="gitter_icon"): lonmin = -23.5 - lonmax = 45.0 - latmin = 29.5 - latmax = 70.5 - xsize=np.size(np.arange(lonmin, lonmax+d_grad, d_grad)) - ysize=np.size(np.arange(latmin, latmax+d_grad, d_grad)) - test_string='#\n# gridID 1\n#\ngridtype = lonlat\ngridsize = %d\nxname = lon\nxlongname = longitude\nxunits = degrees_east\nyname = lat\nylongname = latitude\nyunits = degrees_north\nxsize = %d\nysize = %d\nxfirst = -23.5\nxinc = %.1f\nyfirst = 29.5\nyinc = %.1f' %(xsize*ysize, xsize, ysize, d_grad, d_grad) - fd = open(filename, "x") + lonmax = 45.0 + latmin = 29.5 + latmax = 70.5 + xsize = np.size(np.arange(lonmin, lonmax+d_grad, d_grad)) + ysize = np.size(np.arange(latmin, latmax+d_grad, d_grad)) + test_string = (f"#\n# gridID 1\n#\ngridtype = lonlat\ngridsize = {xsize*ysize}\nxname = lon\nxlongname = longitude" + f"\nxunits = degrees_east\nyname = lat\nylongname = latitude\nyunits = degrees_north" + f"\nxsize = {xsize}\nysize = {ysize}\nxfirst = -23.5\nxinc = {d_grad}\n" + f"yfirst = 29.5\nyinc = {d_grad}" + ) + fd = open(filename, "x") fd.write(test_string) fd.close() + creategridinfo() creategridinfo(d_grad=0.12, filename="gitter_icon1.2") creategridinfo(d_grad=0.2, filename="gitter_icon2") diff --git a/src/main.py b/src/main.py index 95607b9..c27372a 100644 --- a/src/main.py +++ b/src/main.py @@ -12,73 +12,22 @@ import plotlib import modelinfolib as model model = model.icon_nest -# --------------------------------------------------------------------------------------------------------------------- - -def main(): - config = ut.load_yaml('config.yml') - - # command line arguments - parser = argparse.ArgumentParser(description="Hodograph Maps") - # Add positional argument - parser.add_argument('mode', type=str, - help='Mode: Test, Basic, Sounding, or Nixon') - - # Add optional argument - parser.add_argument('-f', '--field', type=str, - help='Which background field: CAPE ML, CAPE CON or LPI ') - - parser.add_argument('-d', '--date', type=str, - help='Date Format YYYY-MM-DD') - - parser.add_argument('-fp', '--fp', type=str, - help='Leadtime or forecast periode') - - parser.add_argument('-r', '--run', type=str, - help='Run (0,6,12,18,.. .etc)') - - # Parse the command line arguments - args = parser.parse_args() - - if args.field == None: - args.field = "CAPE ML" - - if (args.field != "CAPE ML") and (args.field != "CAPE CON") and (args.field != "LPI") and (args.field != "WMAXSHEAR"): - print("Unknown field") - exit(-1) - - if args.date == None: - rundate = datetime.strptime(config["default_date"], "%Y-%m-%d") - else: - rundate = datetime.strptime(args.date, "%Y-%m-%d") - - if args.fp == None: - fp = config["fp"] - else: - fp = args.fp - - if args.run == None: - run = config["run"] - else: - run = args.run - - # replace space with underscores - fieldname = args.field.replace(" ", "_") +# --------------------------------------------------------------------------------------------------------------------- - print(f"\nDate: {rundate}\n Arguments: {args} \nConfig-File: {config}\n\n") - #ut.download_nwp(fieldname, datum="20240227", run="00", fp=0, store_path="./") +def run(program_mode, fieldname, rundate, model_run, fp): + config = ut.load_yaml('config.yml') - - if args.mode == "Test": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, run, fp, path="./iconnest/") + 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, run, titel='CAPE') + plotlib.test_plot(cape_fld, lats, lons, fp, model_run, titel='CAPE') - elif args.mode == "Sounding": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, run, fp, path="./iconnest/") + elif program_mode == "Sounding": + cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./iconnest/") - steps=config["steps"] + steps = config["steps"] nlvl = int(model.getnlev()/steps) t_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) t_fld.fill(np.nan) @@ -95,20 +44,20 @@ def main(): 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, run, fp, path="./iconnest/") - q_fld[lvl_idx,:,:] = ut.open_gribfile_multi("QV", level, rundate, run, fp, path="./iconnest/") - p_fld[lvl_idx,:,:] = ut.open_gribfile_multi("P", level, rundate, run, fp, path="./iconnest/") + 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/") 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, run, titel='CAPE') - elif args.mode == "Basic": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, run, fp, path="./iconnest/") + 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') + elif program_mode == "Basic": + cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./iconnest/") - steps=config["steps"] + steps = config["steps"] nlvl = int(model.getnlev()/steps) u_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) u_fld.fill(np.nan) @@ -123,20 +72,22 @@ def main(): 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, run, fp, path="./iconnest/") - v_fld[lvl_idx,:,:] = ut.open_gribfile_multi("V", level, rundate, run, fp, path="./iconnest/") + 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/") 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, run, titel='CAPE', threshold=config["threshold"]) - plotlib.basic_plot_custarea (cape_fld, u_fld, v_fld, lats, lons, fp, run, titel='CAPE', threshold=config["threshold"]) - elif args.mode == "Nixon": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, run, fp, path="./iconnest/") - - steps=config["steps"] + + 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"]) + elif program_mode == "Nixon": + cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./iconnest/") + + steps = config["steps"] nlvl = int(model.getnlev()/steps) u_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) u_fld.fill(np.nan) @@ -155,26 +106,90 @@ def main(): 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, run, fp, path="./iconnest/") - v_fld[lvl_idx,:,:] = ut.open_gribfile_multi("V", level, rundate, run, fp, path="./iconnest/") - h_fld[lvl_idx,:,:] = ut.open_gribfile_multi("H", level, rundate, run, fp, path="./iconnest/") - p_fld[lvl_idx,:,:] = ut.open_gribfile_multi("P", level, rundate, run, fp, path="./iconnest/") + 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/") lvl_idx += 1 if lvl_idx >= nlvl: break - print(np.nanmean(p_fld, axis=(1,2))) + print(np.nanmean(p_fld, axis=(1, 2))) - du = np.subtract(u_fld[30,:,:], u_fld[0,:,:]) - dv = np.subtract(v_fld[30,:,:], v_fld[0,:,:]) + 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, run, imfmt="png") + plotlib.nixon_proj(cape_fld, dls_fld, u_fld, v_fld, p_fld, h_fld, lats, lons, fp, model_run, imfmt="png") else: print("Wrong command line argument") exit(-1) +# --------------------------------------------------------------------------------------------------------------------- + + +def main(): + config = ut.load_yaml('config.yml') + + # command line arguments + parser = argparse.ArgumentParser(description="Hodograph Maps") + # Add positional argument + parser.add_argument('mode', type=str, + help='Mode: Test, Basic, Sounding, or Nixon') + + # Add optional argument + parser.add_argument('-f', '--field', type=str, + help='Which background field: CAPE ML, CAPE CON or LPI ') + + parser.add_argument('-d', '--date', type=str, + help='Date Format YYYY-MM-DD') + + parser.add_argument('-fp', '--fp', type=str, + help='Leadtime or forecast periode') + + parser.add_argument('-r', '--run', type=str, + help='Run (0,6,12,18,.. .etc)') + + # Parse the command line arguments + args = parser.parse_args() + + if args.field is None: + args.field = "CAPE ML" + + if (args.field != "CAPE ML") and (args.field != "CAPE CON") and (args.field != "LPI") and (args.field != "WMAXSHEAR"): + print("Unknown field") + exit(-1) + + if args.date is None: + rundate = datetime.strptime(config["default_date"], "%Y-%m-%d") + else: + rundate = datetime.strptime(args.date, "%Y-%m-%d") + + if args.fp is None: + fp = config["fp"] + else: + fp = args.fp + + if args.run is None: + model_run = config["run"] + else: + model_run = args.run + + if args.mode != "Test" or args.mode != "Sounding" or args.mode != "Basic" or args.mode != "Nixon": + print("Unknown Mode. Exit program.") + exit(0) + else: + program_mode = args.mode + + # 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) + +# --------------------------------------------------------------------------------------------------------------------- + if __name__ == "__main__": main() - diff --git a/src/meteolib.py b/src/meteolib.py index a8ca319..83f8778 100644 --- a/src/meteolib.py +++ b/src/meteolib.py @@ -1,7 +1,7 @@ #!/usr/bin/python3 ################################ # -# Meteorology library +# Meteorology library # ################################ @@ -15,11 +15,12 @@ GRAV = 9.80665 # Gravity TOL = 1e-10 # Floating Point Tolerance -def uv2spddir(u,v): - direction = np.rad2deg(np.arctan2(-u,v)) +def uv2spddir(u, v): - if type(direction) == np.ndarray: + direction = np.rad2deg(np.arctan2(-u, v)) + + if isinstance(direction, np.ndarray): direction[np.where(direction < 0)] += 360 else: if direction < 0: @@ -39,6 +40,7 @@ def q_to_mixrat(q): """ return q/(1.0-q) + def temp_at_mixrat(w, p): """ Returns the temperature in K of air at the given mixing ratio in g/kg and pressure in hPa @@ -55,11 +57,15 @@ def temp_at_mixrat(w, p): """ # Constants Used - c1 = 0.0498646455 ; c2 = 2.4082965 ; c3 = 7.07475 - c4 = 38.9114 ; c5 = 0.0915 ; c6 = 1.2035 + c1 = 0.0498646455 + c2 = 2.4082965 + c3 = 7.07475 + c4 = 38.9114 + c5 = 0.0915 + c6 = 1.2035 x = np.log10(w * p / (622. + w)) - x = (np.power(10.,((c1 * x) + c2)) - c3 + (c4 * np.power((np.power(10,(c5 * x)) - c6),2))) + x = (np.power(10., ((c1 * x) + c2)) - c3 + (c4 * np.power((np.power(10, (c5 * x)) - c6), 2))) return x @@ -92,12 +98,12 @@ def non_parcel_bunkers_motion_experimental(u, v, ps, i_500m, i_5km, i_6km): """ - d=7.5 + d = 7.5 # sfc-500m Mean Wind mnu500m, mnv500m = mean_wind(u[:i_500m], v[:i_500m], ps[:i_500m]) - ## 5.5km-6.0km Mean Wind + # 5.5km-6.0km Mean Wind mnu5500m_6000m, mnv5500m_6000m = mean_wind(u[i_5km:i_6km], v[i_5km:i_6km], ps[i_5km:i_6km]) # shear vector of the two mean winds @@ -105,7 +111,7 @@ def non_parcel_bunkers_motion_experimental(u, v, ps, i_500m, i_5km, i_6km): shrv = mnv5500m_6000m - mnv500m # SFC-6km Mean Wind - mnu6, mnv6 = mean_wind(u[:i_6km], v[:i_6km], ps[:i_6km]) + mnu6, mnv6 = mean_wind(u[:i_6km], v[:i_6km], ps[:i_6km]) # Bunkers Right Motion tmp = d / np.sqrt(shru*shru + shrv*shrv) @@ -114,4 +120,4 @@ def non_parcel_bunkers_motion_experimental(u, v, ps, i_500m, i_5km, i_6km): lstu = mnu6 - (tmp * shrv) lstv = mnv6 + (tmp * shru) - return rstu, rstv, lstu, lstv, mnu6, mnv6 \ No newline at end of file + return rstu, rstv, lstu, lstv, mnu6, mnv6 diff --git a/src/modelinfolib.py b/src/modelinfolib.py index f34276a..c81012d 100644 --- a/src/modelinfolib.py +++ b/src/modelinfolib.py @@ -1,6 +1,6 @@ #!/usr/bin/python3 """ -## COSMO-REA2 +# COSMO-REA2 points = 564720 #(724x780) nlon = 724 nlat = 780 @@ -16,7 +16,7 @@ #noth_lat= -## COSMO-RE6 +# COSMO-RE6 points = 698752 #(848x824) nlon = 848 nlat = 824 @@ -32,7 +32,7 @@ #noth_lat= -## ICON Nest +# ICON Nest points = 904689 #(1377x657) nlon = 1377 nlat = 657 @@ -45,7 +45,7 @@ d_grad = 0.0625 -## COSMO D2 +# COSMO D2 points = 466116 #(651x716) nlon = 651 nlat = 716 @@ -54,13 +54,13 @@ lonmax = 5.5 latmin = -6.3 latmax = 8.0 -##gradients +#gradients d_grad = 0.02 #nortpole north_lon=-170 noth_lat=40 -##ERAINTERIM +# ERAINTERIM points=115680 #(480x241) nlon=480 nlat=241 @@ -119,7 +119,8 @@ def getlevtyp(self): def getd_grad(self): return self.d_grad + # Example usage: -icon_nest = MODELIFNO(1377, 657, 0.0625, 54, "model") # lowest 74 and we download till 20 +icon_nest = MODELIFNO(1377, 657, 0.0625, 54, "model") # lowest 74 and we download till 20 cosmo_d2 = MODELIFNO(651, 716, 0.02, 65, "model") ifs = MODELIFNO(450, 900, 0.4, 10, "pres") diff --git a/src/plotlib.py b/src/plotlib.py index 7e96fc5..60351a3 100644 --- a/src/plotlib.py +++ b/src/plotlib.py @@ -20,7 +20,7 @@ # own moduls import utilitylib as ut import meteolib as met -from skewT import * +from skewT import SkewXAxes # --------------------------------------------------------------------------------------------------------------------- # create plot class @@ -32,6 +32,7 @@ fontsize = config["fontsize"] titlesize = config["titlesize"] + def eu_merc(hour, start, 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]) @@ -98,6 +99,7 @@ def alps(hour, start, projection=crs.EuroPP(), lon1=5.8, lon2=17.8, lat1=45.23, plt.annotate(string2, xy=(0, 1), xycoords='axes fraction', fontsize=fontsize) return fig, ax + def two_plots(projection=crs.EuroPP(), lon1=3.56, lon2=16.5, lat1=46.2, lat2=55.6, fac=3): fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(5*fac, 3*fac), subplot_kw=dict(projection=projection)) fig.subplots_adjust(left=0.02, right=0.92, top=0.95, bottom=0.05, wspace=0.14) @@ -110,18 +112,18 @@ def two_plots(projection=crs.EuroPP(), lon1=3.56, lon2=16.5, lat1=46.2, lat2=55. plt.annotate("ICON Nest (DWD)", xy=(0.6, -0.03), xycoords='axes fraction', fontsize=fontsize) return fig, ax1, ax2 -# --------------------------------------------------------------------------------------------------------------------- +# --------------------------------------------------------------------------------------------------------------------- # create colormap for CAPE field clevs = np.array([50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000]) cmap = LinearSegmentedColormap.from_list("", ["green", "yellow", "orange", "red", "darkred", "darkmagenta"]) -cmap2= LinearSegmentedColormap.from_list("", ["gold", "orange", "darkorange", "red", "darkred", "darkmagenta"]) +cmap2 = LinearSegmentedColormap.from_list("", ["gold", "orange", "darkorange", "red", "darkred", "darkmagenta"]) # --------------------------------------------------------------------------------------------------------------------- -def test_plot (cape_fld, lats, lons, hour, run, titel='CAPE'): +def test_plot(cape_fld, lats, lons, hour, run, titel='CAPE'): """ Parameters: ------------ @@ -140,8 +142,8 @@ def test_plot (cape_fld, lats, lons, hour, run, titel='CAPE'): plt.title(titel, fontsize=titlesize) wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), - cmap=cmap, extend = 'max', alpha=0.4, antialiased=True) - + cmap=cmap, extend='max', alpha=0.4, antialiased=True) + cax = fig.add_axes([0.27, 0.05, 0.35, 0.05]) fig.colorbar(wx, cax=cax, orientation='horizontal') @@ -151,6 +153,7 @@ def test_plot (cape_fld, lats, lons, hour, run, titel='CAPE'): # --------------------------------------------------------------------------------------------------------------------- + def soundingpoint(point, temp, q, pres, ax, width=0.1, proj='skewT', smooth=False): """ Parameters: @@ -183,26 +186,27 @@ def soundingpoint(point, temp, q, pres, ax, width=0.1, proj='skewT', smooth=Fals # ax2.set_ylim(0, ) # 10 ms circle - #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)+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) # calculate dewpoint pres /= 100 - dewpoint = met.temp_at_mixrat( met.q_to_mixrat(q)*1000.0, pres ) + dewpoint = met.temp_at_mixrat(met.q_to_mixrat(q)*1000.0, pres) # smoothing if smooth is True: temp[1:] = (temp[1:] + temp[:-1])/2 dewpoint[1:] = (dewpoint[1:] + dewpoint[:-1])/2 - + ax2.semilogy(temp[:40:1]-met.ZEROCNK, pres[:40:1], 'b-', lw=1.5) ax2.semilogy(dewpoint[:40:1]-met.ZEROCNK, pres[:40:1], 'g-', lw=1.5) - ax2.vlines(x=0, ymin=300, ymax=1000, colors='purple', ls='--', lw=0.6) # freezing level - ax2.set_ylim(1050,300) - #ax2.invert_xaxis() + ax2.vlines(x=0, ymin=300, ymax=1000, colors='purple', ls='--', lw=0.6) # freezing level + ax2.set_ylim(1050, 300) + # ax2.invert_xaxis() # --------------------------------------------------------------------------------------------------------------------- + def sounding_plot(cape_fld, temp, q_fld, p_fld, lats, lons, hour, start, titel='CAPE', imfmt="png"): """ Parameters: @@ -222,19 +226,16 @@ def sounding_plot(cape_fld, temp, q_fld, p_fld, lats, lons, hour, start, titel=' fig, ax = ce_states(hour, start, projection=crs.PlateCarree()) plt.title(titel, fontsize=titlesize) - wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), cmap=cmap2, + wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), cmap=cmap2, extend='max', alpha=0.4, antialiased=True) - - - for i in range(275, 415, 10): for j in range(420, 670, 15): soundingpoint((lons[i, j], lats[i, j]), - 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() + 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() cax = fig.add_axes([0.27, 0.05, 0.35, 0.05]) fig.colorbar(wx, cax=cax, orientation='horizontal') @@ -248,6 +249,7 @@ def sounding_plot(cape_fld, temp, q_fld, p_fld, lats, lons, hour, start, titel=' # --------------------------------------------------------------------------------------------------------------------- + def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): """ Parameters: @@ -262,7 +264,6 @@ def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): -------- None """ - # this takes us from the data coordinates to the display coordinates. test = ax.transData.transform(point) @@ -297,7 +298,7 @@ def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): 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), - np.zeros(100)+30, '-k', alpha=.3, lw=0.8) + np.zeros(100)+30, '-k', alpha=.3, lw=0.8) ax2.plot(wdir[:10:1], spd[:10:1], 'r-', lw=1.5) ax2.plot(wdir[9:21:2], spd[9:21:2], 'g-', lw=1.5) @@ -326,7 +327,7 @@ def basic_plot(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold= fig, ax = ce_states(hour, start, projection=crs.PlateCarree()) plt.title(titel, fontsize=titlesize) - wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), cmap=cmap, + 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): @@ -351,9 +352,9 @@ def basic_plot(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold= name = f"./images/hodographmap_ce_{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(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold=10., imfmt="png"): lon1 = config["customize"]["lon1"] lon2 = config["customize"]["lon2"] lat1 = config["customize"]["lat1"] @@ -361,7 +362,8 @@ def basic_plot_custarea(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', t fig, ax = customize_area(hour, start, projection=crs.PlateCarree(), lon1=lon1, lon2=lon2, lat1=lat1, lat2=lat2) plt.title(titel, fontsize=titlesize) - wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), cmap=cmap, extend = 'max', alpha=0.4, antialiased=True) + 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(340, 400, 5): for j in range(555, 665, 5): @@ -395,7 +397,7 @@ def nixon_hodograph(point, u, v, p, height, ax, width=0.1, clim=40, proj='polar' rstu, rstv : storm motion vector for right mover """ - i=0 + i = 0 while height[i] < 500: i += 1 i_500m = deepcopy(i) @@ -409,8 +411,8 @@ def nixon_hodograph(point, u, v, p, height, ax, width=0.1, clim=40, proj='polar' rstu, rstv, lstu, lstv, mwu6, mwv6 = met.non_parcel_bunkers_motion_experimental(u, v, p, i_500m, i_5km, i_6km) - u-=rstu - v-=rstv + u -= rstu + v -= rstv # plot test = ax.transData.transform(point) @@ -442,7 +444,7 @@ def nixon_hodograph(point, u, v, p, height, ax, width=0.1, clim=40, proj='polar' ax2.plot(wdir[9:21:2], spd[9:21:2], 'g-', lw=1.5) ax2.plot(wdir[19:-20:2], spd[19:-20:2], 'b-', lw=1.5) ax2.scatter(0, 0, c="k", s=2, marker='x', alpha=0.75) - + theta, mag = met.uv2spddir(rstu, rstv) ax2.arrow(theta, 0, 0, mag, head_width=0.1, head_length=0.1) @@ -454,9 +456,8 @@ def nixon_proj(cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt= dls_fld : 2D deep layer shear field or brn shear ... u, v : wind components high : model level high - or + or rstu, rstv : storm motion vector - background filed is cape ... only hodographs with more the 10 m/s dls """ @@ -468,7 +469,7 @@ def nixon_proj(cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt= 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]) - cmap = LinearSegmentedColormap.from_list("", ["green", "yellow", "orange", "red", "darkred", "darkmagenta"]) + cmap = LinearSegmentedColormap.from_list("", ["green", "yellow", "orange", "red", "darkred", "darkmagenta"]) wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), cmap=cmap, extend='max', alpha=0.4) @@ -478,10 +479,10 @@ def nixon_proj(cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt= # cleves = np.array([500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000]) # cs = plt.contour(lons, lats, wx_fld[:, :], levels=cleves, transform=crs.PlateCarree(), colors='k', linewidths=0.8) # plt.clabel(cs, np.array([500, 1000, 2000, 3000]), fontsize=7, inline=1, fmt='%.f') # contour labels - + for i in range(280, 410, 10): for j in range(420, 670, 15): - if np.mean(dls_fld[i-1:i+1, j-1:j+1]) > 10.0 and np.mean(cape_fld[i-1:i+1, j-1:j+1]) > 10.0:#m/s + if np.mean(dls_fld[i-1:i+1, j-1:j+1]) > 10.0 and np.mean(cape_fld[i-1:i+1, j-1:j+1]) > 10.0: # m/s nixon_hodograph((lons[i, j], lats[i, j]), np.mean(u[::-1, i-1:i+1, j-1:j+1], axis=(1, 2)), np.mean(v[::-1, i-1:i+1, j-1:j+1], axis=(1, 2)), @@ -491,7 +492,7 @@ def nixon_proj(cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt= cax = fig.add_axes([0.27, 0.05, 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) @@ -499,4 +500,4 @@ def nixon_proj(cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt= name = f"./images/nixon_ce_{hour}.{imfmt}" plt.savefig(name) - plt.close() \ No newline at end of file + plt.close() diff --git a/src/skewT.py b/src/skewT.py index 8ff84fb..eebb494 100644 --- a/src/skewT.py +++ b/src/skewT.py @@ -2,7 +2,7 @@ import matplotlib.transforms as transforms import matplotlib.axis as maxis import matplotlib.spines as mspines -from matplotlib.projections import register_projection +# from matplotlib.projections import register_projection # The sole purpose of this class is to look at the upper, lower, or total diff --git a/src/utilitylib.py b/src/utilitylib.py index 46beed0..ee4dae0 100644 --- a/src/utilitylib.py +++ b/src/utilitylib.py @@ -9,7 +9,7 @@ # --------------------------------------------------------------------------------------------------------------------- -def datum (h, start): +def datum(h, start): if not isinstance(h, int): h = int(h) @@ -23,7 +23,7 @@ def datum (h, start): 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") - + return string1, string2 # --------------------------------------------------------------------------------------------------------------------- @@ -45,15 +45,14 @@ def load_yaml(yaml_file, yaml_path='./'): def download_nwp(fieldname, datum="20240227", run="00", fp=0, store_path="./"): - - opendataserver="https://opendata.dwd.de/weather/nwp/icon-eu/grib" - nwp_name=f"icon-eu_europe_regular-lat-lon_single-level_{datum}{run}_{fp:03d}_{fieldname.upper()}.grib2.bz2" - url_link=f"{opendataserver}/{run}/{fieldname.lower()}/{nwp_name}" + opendataserver = "https://opendata.dwd.de/weather/nwp/icon-eu/grib" + nwp_name = f"icon-eu_europe_regular-lat-lon_single-level_{datum}{run}_{fp:03d}_{fieldname.upper()}.grib2.bz2" + url_link = f"{opendataserver}/{run}/{fieldname.lower()}/{nwp_name}" response = requests.get(url_link) with open(f"{store_path}/test.grib2.bz2", 'wb') as f: f.write(response.content) - + print("Download complete.") # --------------------------------------------------------------------------------------------------------------------- @@ -80,6 +79,7 @@ def open_gribfile_single(fieldname, datetime_obj, run, fp, path="./iconnest/"): print("Maximum:", np.nanmax(data)) return data, lats, lons + def open_gribfile_multi(fieldname, lvl, datetime_obj, run, fp, path="./iconnest/"): date_string = datetime_obj.strftime("%Y%m%d") nwp_modellevel = "icon-eu_europe_regular-lat-lon_model-level"