diff --git a/.gitignore b/.gitignore index 18d9c32..96cb0ae 100644 --- a/.gitignore +++ b/.gitignore @@ -2,10 +2,6 @@ .ipynb_checkpoints/ */.ipynb_checkpoints/ -*/*/.ipynb_checkpoints/ - -*/*.csv -*/*/*.csv */*.png */*/*.png @@ -13,11 +9,12 @@ */*.svg */*/*.svg -*/*.jpg -*/*/*.jpg - *.log */*.log +*.pyc +*/*.pyc +*/*/*.pyc + src/invariant src/iconnest diff --git a/README.md b/README.md index 1560323..b4a174a 100644 --- a/README.md +++ b/README.md @@ -11,14 +11,26 @@ conda activate HodographMaps ``` -# Example Image +## Example Image -![Example](images/test_ce_9.png) +![Example](images/example.png) +![Example](images/hodographmap_ce_9.png) +## How to run + ```bash cd src bash download_script.sh conda activate HodographMaps ``` + +## Cartopy? + +https://github.com/mammatus95/cartopy +https://scitools.org.uk/cartopy/docs/latest/# + +## Datasource +ICON Nest: https://opendata.dwd.de/weather/nwp/icon-eu/ +IFS: https://www.ecmwf.int/en/forecasts/datasets/open-data diff --git a/images/example.png b/images/example.png new file mode 100644 index 0000000..8f42aa7 Binary files /dev/null and b/images/example.png differ diff --git a/images/hodographmap_ce_9.png b/images/hodographmap_ce_9.png new file mode 100644 index 0000000..92bcac0 Binary files /dev/null and b/images/hodographmap_ce_9.png differ diff --git a/images/test_ce_9.png b/images/test_ce_9.png deleted file mode 100644 index 5ba961d..0000000 Binary files a/images/test_ce_9.png and /dev/null differ diff --git a/src/__pycache__/plotlib.cpython-311.pyc b/src/__pycache__/plotlib.cpython-311.pyc deleted file mode 100644 index c9040bf..0000000 Binary files a/src/__pycache__/plotlib.cpython-311.pyc and /dev/null differ diff --git a/src/__pycache__/utilitylib.cpython-311.pyc b/src/__pycache__/utilitylib.cpython-311.pyc deleted file mode 100644 index e682a94..0000000 Binary files a/src/__pycache__/utilitylib.cpython-311.pyc and /dev/null differ diff --git a/src/config.yml b/src/config.yml new file mode 100644 index 0000000..b3414f2 --- /dev/null +++ b/src/config.yml @@ -0,0 +1,12 @@ +run: 0 +fp: 9 +default_date: "2024-02-27" + +levels: [60,10] +steps: 1 # level steps + +threshold: 0 + +fontsize: 12 +titlesize: 18 + diff --git a/src/download_script.sh b/src/download_script.sh index 8dcd3fa..af0e4b0 100755 --- a/src/download_script.sh +++ b/src/download_script.sh @@ -63,7 +63,7 @@ do cd ${path_icon} T=$(printf "%03d" "$X") #single level - for N in CAPE_ML + for N in CAPE_ML #CAPE_CON do typeset -l nvar nvar=${N} @@ -83,7 +83,7 @@ do done #rm -f - + python3 main.py Basic --fp ${X} --run ${R} --date $(date) echo "done with ${T}" done diff --git a/src/main.py b/src/main.py index b9c25b7..c537f24 100644 --- a/src/main.py +++ b/src/main.py @@ -1,15 +1,21 @@ #!/usr/bin/python3 +from datetime import datetime + + import argparse +import numpy as np + # own moduls import utilitylib as ut import plotlib - +import modelinfolib as model +model = model.icon_nest # --------------------------------------------------------------------------------------------------------------------- def main(): - # config = ul.load_yaml('config.yml') + config = ut.load_yaml('config.yml') # command line arguments parser = argparse.ArgumentParser(description="Hodograph Maps") @@ -20,6 +26,15 @@ def main(): # 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() @@ -31,20 +46,66 @@ def main(): 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(args) + print(f"\nDate: {rundate}\n Arguments: {args} \nConfig-File: {config}\n\n") #ut.download_nwp(fieldname, datum="20240227", run="00", fp=0, store_path="./") - cape_fld, lats, lons = ut.open_gribfile(fieldname, path="./iconnest/") - print(lats.shape) + if args.mode == "Test": - plotlib.test_plot (cape_fld, lats, lons, 9, 0, titel='CAPE') + cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, 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') + else: + cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, 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) + v_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + v_fld.fill(np.nan) + + if config["levels"][0] > config["levels"][1]: + steps *= -1 + elif config["levels"][0] == config["levels"][1]: + print("Wrong levels in config.yml!") + exit(0) + + 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/") + + 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, 9, 0, titel='CAPE', threshold=config["threshold"]) + + if __name__ == "__main__": main() diff --git a/src/modelinfolib.py b/src/modelinfolib.py new file mode 100644 index 0000000..02e433e --- /dev/null +++ b/src/modelinfolib.py @@ -0,0 +1,116 @@ +#!/usr/bin/python3 +""" +## COSMO-REA2 +points = 564720 #(724x780) +nlon = 724 +nlat = 780 +nlev = 50 +lonmin = -7.5 +lonmax = 5.514 +latmin = -6.0 +latmax = 8.022 +##gradients +d_grad = 0.018 +#nortpole +#north_lon= +#noth_lat= + + +## COSMO-RE6 +points = 698752 #(848x824) +nlon = 848 +nlat = 824 +nlev = 40 +lonmin = -28.403 +lonmax = 18.182 +latmin = -23.403 +latmax = 21.862 +##gradients +d_grad = 0.055 +#nortpole +#north_lon= +#noth_lat= + + +## ICON Nest +points = 904689 #(1377x657) +nlon = 1377 +nlat = 657 +nlev = 51 +lonmin = -23.5 +lonmax = 45.0 +latmin = 29.5 +latmax = 70.5 +##gradients +d_grad = 0.0625 + + +## COSMO D2 +points = 466116 #(651x716) +nlon = 651 +nlat = 716 +nlev = 65 +lonmin = -7.5 +lonmax = 5.5 +latmin = -6.3 +latmax = 8.0 +##gradients +d_grad = 0.02 +#nortpole +north_lon=-170 +noth_lat=40 + +##ERAINTERIM +points=115680 #(480x241) +nlon=480 +nlat=241 +nlev=23 +lonmin = -180.0 +lonmax = 179.25 +latmin = -90 +latmax = 90 +d_grad = 0.75 +""" + +class MODELIFNO: + def __init__(self, nlon, nlat, nlev, d_grad): + self.points = nlon*nlat + self.nlon = nlon + self.nlat = nlat + self.nlev = nlev + self.d_grad = d_grad + + 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) + + def __str__(self): + return ( + f"Model Information:\nPoints: {self.points}\n" + f"Number of Longitudes: {self.nlon}\nNumber of Latitudes: {self.nlat}\n" + f"Number of Levels: {self.nlev}\nGradient: {self.d_grad}" + ) + + def getpoints(self): + return self.points + + def getnlon(self): + return self.nlon + + def getnlat(self): + return self.nlat + + def getnlev(self): + return self.nlev + + def getd_grad(self): + return self.d_grad + +# Example usage: +icon_nest = MODELIFNO(1377, 657, 51, -23.5, 45.0, 29.5, 70.5, 0.0625) +cosmo_d2 = MODELIFNO(651, 716, 65, 0.02) +ifs = MODELIFNO(450, 900, 10, 0.4) diff --git a/src/plotlib.py b/src/plotlib.py index 57ff2e1..02f4ebe 100644 --- a/src/plotlib.py +++ b/src/plotlib.py @@ -8,9 +8,8 @@ import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt -from matplotlib.colors import ListedColormap, LinearSegmentedColormap +from matplotlib.colors import ListedColormap, BoundaryNorm, LinearSegmentedColormap from matplotlib.cm import get_cmap -from matplotlib.colors import BoundaryNorm from matplotlib.ticker import MaxNLocator from mpl_toolkits.axes_grid1 import make_axes_locatable @@ -30,11 +29,15 @@ # own moduls import utilitylib as ut import kinematiclib as kin + # --------------------------------------------------------------------------------------------------------------------- # create plot class # ===================================================================================================================== # https://scitools.org.uk/cartopy/docs/v0.16/crs/projections.html#cartopy-projections +config = ut.load_yaml('config.yml') +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)) @@ -47,8 +50,8 @@ def eu_merc(hour, start, projection=crs.Mercator(), factor=3): 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=10) - plt.annotate(string2, xy=(0, 1), 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) return fig, ax @@ -60,8 +63,8 @@ def eu_states(hour, start, projection=crs.EuroPP()): 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=10) - plt.annotate(string2, xy=(0, 1), 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) return fig, ax def ce_states(hour, start, projection=crs.EuroPP(), lon1=1.56, lon2=18.5, lat1=45.1, lat2=56.6): @@ -72,8 +75,8 @@ def ce_states(hour, start, projection=crs.EuroPP(), lon1=1.56, lon2=18.5, lat1=4 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), xycoords='axes fraction', fontsize=10) - plt.annotate(string2, xy=(0, 1), 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) return fig, ax def east_de(hour, start, projection=crs.EuroPP(), lon1=10.7, lon2=18, lat1=49.8, lat2=54.8): @@ -85,8 +88,8 @@ def east_de(hour, start, projection=crs.EuroPP(), lon1=10.7, lon2=18, lat1=49.8, 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), xycoords='axes fraction', fontsize=10) - plt.annotate(string2, xy=(0, 1), xycoords='axes fraction', fontsize=10) + plt.annotate(string1, xy=(0.835, 1), xycoords='axes fraction', fontsize=fontsize) + plt.annotate(string2, xy=(0, 1), 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): @@ -96,8 +99,8 @@ def alps(hour, start, projection=crs.EuroPP(), lon1=5.8, lon2=17.8, lat1=45.23, 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=10) - plt.annotate(string2, xy=(0, 1), 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) return fig, ax def two_plots(projection=crs.EuroPP(), lon1=3.56, lon2=16.5, lat1=46.2, lat2=55.6, fac=3): @@ -109,12 +112,20 @@ def two_plots(projection=crs.EuroPP(), lon1=3.56, lon2=16.5, lat1=46.2, lat2=55. ax2.set_extent([lon1, lon2, lat1, lat2]) ax2.coastlines('10m', linewidth=1.2) ax2.add_feature(states_provinces, edgecolor='black') - plt.annotate("ICON Nest (DWD)", xy=(0.6, -0.03), xycoords='axes fraction', fontsize=10) + 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"]) + +# --------------------------------------------------------------------------------------------------------------------- + + def test_plot (cape_fld, lats, lons, hour, run, titel='CAPE'): """ Parameters: @@ -129,13 +140,9 @@ def test_plot (cape_fld, lats, lons, hour, run, titel='CAPE'): -------- None """ - - # create colormap for CAPE field - clevs = np.array([1, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000]) - cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", ["green", "yellow", "orange", "red", "darkred", "darkmagenta"]) fig, ax = ce_states(hour, run, projection=crs.PlateCarree()) - plt.title(titel, fontsize=15) + 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) @@ -146,72 +153,67 @@ def test_plot (cape_fld, lats, lons, hour, run, titel='CAPE'): plt.savefig(name) plt.close() - # --------------------------------------------------------------------------------------------------------------------- -def hodopoint(point, u, v, ax, width=0.1, clim=40, projection='polar'): +def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): """ Parameters: ------------ - width : width of added axes for the hodo as fig coordinates (examples: 0.5 are half of the image) - clim : max. wind speed magnitude of the hodo - point : tuple of data coordinates (10, 55) - u, v : wind components - projection : 'polar' + width : width of added axes for the hodo as fig coordinates (examples: 0.5 are half of the image) + clim : max. wind speed magnitude of the hodo + point : tuple of data coordinates (10, 55) + u, v : wind components + proj : 'polar' Returns: -------- None """ - #print(point) - #proj_cart = crs.PlateCarree() - #point=(0.5, 0.5) - # convert from Axes coordinates to display coordinates - #p_a_disp = ax.transAxes.transform(point) - - # convert from display coordinates to data coordinates - #p_a_data = ax.transData.inverted().transform(p_a_disp) - #this takes us from the data coordinates to the display coordinates. + # this takes us from the data coordinates to the display coordinates. test = ax.transData.transform(point) - #print(test) - #this should take us from the display coordinates to the axes coordinates. + + # this should take us from the display coordinates to the axes coordinates. trans = ax.transAxes.inverted().transform(test) - #print(p_a_data) - # convert from data to cartesian coordinates - #trans = proj_cart.transform_point(src_crs=proj, x=point[0], y=point[1]) - #print(trans) - #trans=(0.5, 0.5) - #create new axe - ax2 = plt.axes([trans[0]-width/2, trans[1]-width/2, width, width], projection=projection) + + # create new axes + ax2 = plt.axes([trans[0]-width/2, trans[1]-width/2, width, width], projection=proj) ax2.get_xaxis().set_visible(False) ax2.get_yaxis().set_visible(False) - #ax2.patch.set_visible(False) ax2.set_frame_on(False) - #ax2.set_xlim(-clim, clim) + + # ax2.set_xlim(-clim, clim) ax2.set_ylim(0, clim) ax2.set_theta_offset(np.pi/2) - #ax2.set_theta_direction(-1) - #10 ms circle + # ax2.set_theta_direction(-1) + + # 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)+30, '-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 = kin.uv2spddir(u, v) - #print("u:", u, "v:", v, " \twdir:", wdir*180./np.pi) - #wdir[1:] = (wdir[1:] + wdir[:-1])/2 - #spd[1:] = (spd[1:] + spd[:-1])/2 + + # smoothing + if smooth == True: + wdir[1:] = (wdir[1:] + wdir[:-1])/2 + spd[1:] = (spd[1:] + spd[:-1])/2 + + # draw part of second cricle if np.max(spd[:-20]) > 28: - kp = np.abs(np.max(wdir[np.where(spd[:-20]>25)])-np.min(wdir[np.where(spd[:-20]>25)]))#/2 - #kp[np.where(kp[:] <= np.pi/8)] = np.pi/8 - 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) + 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) + 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) ax2.plot(wdir[19:-20:2], spd[19:-20:2], 'b-', lw=1.5) ax2.scatter(0, 0, c="k", s=10, marker='x', alpha=0.75) -def basic_plot (cape_fld, u, v, lats, lons, hour, start, titel='CAPE'): + +def basic_plot (cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold=10., imfmt="png"): """ Parameters: ------------ @@ -226,32 +228,20 @@ def basic_plot (cape_fld, u, v, lats, lons, hour, start, titel='CAPE'): -------- None """ - - # 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 = matplotlib.colors.LinearSegmentedColormap.from_list("", ["green", "yellow", "orange", "red", "darkred", "darkmagenta"]) fig, ax = ce_states(hour, start, projection=crs.PlateCarree()) - plt.title(titel, fontsize=15) - - wx = ax.contourf(lons, lats, cape_fld[:, :], levels=clevs, transform=crs.PlateCarree(), cmap=cmap, extend = 'max', alpha=0.4, antialiased=True) - #cb = plt.colorbar(wx, ticks=clevs, shrink=.8) - #cb.set_label(r'$m^2/s^2$') + plt.title(titel, fontsize=titlesize) - #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 + 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(280, 410, 10): - #j=592 + for i in range(275, 415, 10): for j in range(420, 670, 15): - if np.mean(cape_fld[i-1:i+1, j-1:j+1]) > 10.0: #np.mean(ko_fld[i-1:i+1, j-1:j+1]) < 0.0: # - hodopoint((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)), ax, width=0.1, proj=crs.PlateCarree()) - - #divider = make_axes_locatable(ax) - #cax = divider.append_axes("bottom", size="5%", pad=0.01) - #plt.colorbar(wx, cax=cax, orientation='horizontal') - + if np.mean(cape_fld[i-1:i+1, j-1:j+1]) > threshold: + hodopoint((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)), ax, width=0.1) # , proj=crs.PlateCarree() + cax = fig.add_axes([0.27, 0.05, 0.35, 0.05]) fig.colorbar(wx, cax=cax, orientation='horizontal') @@ -261,27 +251,25 @@ def basic_plot (cape_fld, u, v, lats, lons, hour, start, titel='CAPE'): 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) - name = f"hodo_ce_{hour}.png" - plt.savefig(name) - name = f"hodo_ce_{hour}.svg" + + name = f"./images/hodographmap_ce_{hour}.{imfmt}" plt.savefig(name) plt.close() - +def basic_plot_eastde (cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold=10., imfmt="png"): + + fig, ax = east_de(hour, start, projection=crs.PlateCarree()) - plt.title(titel, fontsize=15) + 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) - #cleves = np.array([-20, -2, 20]) - #cs = plt.contour(lons, lats, ko_fld[:, :], levels=cleves, transform=crs.PlateCarree(), colors='k', linewidths=0.8) - #plt.clabel(cs, np.array([-20, -2, 20]), fontsize=7, inline=1, fmt='%.f') - for i in range(340, 400, 5): - #j=592 for j in range(555, 665, 5): - if np.mean(cape_fld[i-1:i+1, j-1:j+1]) > 10.0: #np.mean(ko_fld[i-1:i+1, j-1:j+1]) < -2.0: # - hodopoint((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)), ax, width=0.1, proj=crs.PlateCarree()) + if np.mean(cape_fld[i-1:i+1, j-1:j+1]) > threshold: + hodopoint((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)), ax, width=0.1, proj=crs.PlateCarree()) cax = fig.add_axes([0.27, 0.05, 0.35, 0.05]) fig.colorbar(wx, cax=cax, orientation='horizontal') @@ -292,17 +280,15 @@ def basic_plot (cape_fld, u, v, lats, lons, hour, start, titel='CAPE'): 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) - name = f"hodo_east_{hour}.png" + + name = f"./images/hodographmap_east_{hour}.{imfmt}" plt.savefig(name) - #name = f"hodo_east_{hour}.svg" - #plt.savefig(name) plt.close() # --------------------------------------------------------------------------------------------------------------------- def nixon_hodograph(point, u, v, p, height, ax, width=0.1, clim=40, proj=crs.EuroPP()): - #float(sys.argv[6]), float(sys.argv[7]) """ u, v : horizontal wind components rstu, rstv : storm motion vector for right mover @@ -364,8 +350,7 @@ def nixon_hodograph(point, u, v, p, height, ax, width=0.1, clim=40, proj=crs.Eur ax2.arrow(theta, 0, 0, mag, head_width=0.1, head_length=0.1) -#def nixon_proj (cape_fld, dls_fld, u, v, lat, lon, hour, start): -def nixon_proj (cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start): +def nixon_proj (cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt="png"): """ Nixon projection cape_fld : 2D cape field @@ -378,16 +363,15 @@ def nixon_proj (cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start): background filed is cape ... only hodos with more the 10 m/s dls """ - - + titel = 'CAPE with Hodographs' fig, ax = ce_states(hour, start, projection=crs.PlateCarree()) - plt.title(titel, fontsize=15) + plt.title(titel, 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]) - cmap = matplotlib.colors.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) #cb = plt.colorbar(wx, ticks=clevs, shrink=.8) @@ -401,7 +385,11 @@ def nixon_proj (cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start): #j=592 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 - 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)), 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()) + 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)), + 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]) @@ -412,8 +400,7 @@ def nixon_proj (cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start): 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) - name = f"nixon_ce_{hour}.png" - plt.savefig(name) - name = f"nixon_ce_{hour}.svg" + + name = f"./images/nixon_ce_{hour}.{imfmt}" plt.savefig(name) plt.close() \ No newline at end of file diff --git a/src/utilitylib.py b/src/utilitylib.py index 1e602f4..feeda9f 100644 --- a/src/utilitylib.py +++ b/src/utilitylib.py @@ -37,13 +37,13 @@ def load_yaml(yaml_file, yaml_path='./'): ----------- yaml_file : name of yaml file """ - with open(f"{yaml_path}/{yaml_file}",' r') as yhand: + with open(f"{yaml_path}/{yaml_file}", 'r') as yhand: config_data = yaml.safe_load(yhand) return config_data # --------------------------------------------------------------------------------------------------------------------- # https://opendata.dwd.de/weather/nwp/icon-eu/grib/00/cape_ml/icon-eu_europe_regular-lat-lon_single-level_2024022700_006_CAPE_ML.grib2.bz2 - +# https://opendata.dwd.de/weather/nwp/icon-eu/grib/00/u/icon-eu_europe_regular-lat-lon_model-level_2024022700_000_10_U.grib2.bz2 def download_nwp(fieldname, datum="20240227", run="00", fp=0, store_path="./"): @@ -60,9 +60,12 @@ def download_nwp(fieldname, datum="20240227", run="00", fp=0, store_path="./"): # --------------------------------------------------------------------------------------------------------------------- -def open_gribfile(fieldname, path="./iconnest/"): +def open_gribfile_single(fieldname, datetime_obj, run, fp, path="./iconnest/"): + date_string = datetime_obj.strftime("%Y%m%d") + nwp_singlelevel = "icon-eu_europe_regular-lat-lon_single-level" + # Open the GRIB file - grbs = pygrib.open(f"{path}icon-eu_europe_regular-lat-lon_single-level_2024022700_009_{fieldname.upper()}.grib2") + grbs = pygrib.open(f"{path}{nwp_singlelevel}_{date_string}{run:02d}_{fp:03d}_{fieldname.upper()}.grib2") # Read specific data from the file first_message = grbs[1] @@ -74,10 +77,32 @@ def open_gribfile(fieldname, path="./iconnest/"): lats, lons = first_message.latlons() grbs.close() - print("Data shape:", data.shape) + # print("Data shape:", data.shape) 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" + + # Open the GRIB file + grbs = pygrib.open(f"{path}{nwp_modellevel}_{date_string}{run:02d}_{fp:03d}_{lvl}_{fieldname.upper()}.grib2") + + # Read specific data from the file + first_message = grbs[1] + + data = first_message.values + # Convert missing values to NaNs + data[data == first_message.missingValue] = np.nan + + lats, lons = first_message.latlons() + grbs.close() + + # print("Data shape:", data.shape) + print("Maximum:", np.nanmax(data)) + return data + + # --------------------------------------------------------------------------------------------------------------------- def open_netcdf(fieldname, path="./iconnest/"):