diff --git a/README.md b/README.md index b4a174a..2b0c6c7 100644 --- a/README.md +++ b/README.md @@ -24,13 +24,14 @@ conda activate HodographMaps cd src bash download_script.sh conda activate HodographMaps +python3 main.py Basic ``` ## Cartopy? -https://github.com/mammatus95/cartopy -https://scitools.org.uk/cartopy/docs/latest/# +- 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 +- ICON Nest: https://opendata.dwd.de/weather/nwp/icon-eu/ +- IFS: https://www.ecmwf.int/en/forecasts/datasets/open-data diff --git a/src/config.yml b/src/config.yml index b3414f2..b2a9699 100644 --- a/src/config.yml +++ b/src/config.yml @@ -10,3 +10,11 @@ threshold: 0 fontsize: 12 titlesize: 18 +# customize +customize: # east de + lon1: 10.7 + lon2: 18 + lat1: 49.8 + lat2: 54.8 + +# 10.77, 18.92, 49.80, 55.06 diff --git a/src/main.py b/src/main.py index c537f24..8e81c9c 100644 --- a/src/main.py +++ b/src/main.py @@ -76,7 +76,8 @@ def main(): 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: + + elif args.mode == "Basic": cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, run, fp, path="./iconnest/") steps=config["steps"] @@ -102,9 +103,48 @@ def main(): 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/") - plotlib.basic_plot (cape_fld, u_fld, v_fld, lats, lons, 9, 0, titel='CAPE', threshold=config["threshold"]) + 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) + h_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + h_fld.fill(np.nan) + p_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) + p_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/") + 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/") + + lvl_idx += 1 + if lvl_idx >= nlvl: + break + + 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,:,:]) + 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") + else: + print("Wrong command line argument") + exit(-1) if __name__ == "__main__": diff --git a/src/modelinfolib.py b/src/modelinfolib.py index 02e433e..08053b6 100644 --- a/src/modelinfolib.py +++ b/src/modelinfolib.py @@ -72,7 +72,9 @@ d_grad = 0.75 """ + class MODELIFNO: + def __init__(self, nlon, nlat, nlev, d_grad): self.points = nlon*nlat self.nlon = nlon @@ -80,6 +82,7 @@ def __init__(self, nlon, nlat, nlev, d_grad): 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 @@ -87,6 +90,7 @@ def __init__(self, nlon, nlat, nlev, lonmin, lonmax, latmin, latmax, d_grad): self.nlev = nlev self.d_grad = d_grad print(lonmin, lonmax, latmin, latmax) + """ def __str__(self): return ( @@ -94,13 +98,13 @@ def __str__(self): 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 @@ -110,7 +114,8 @@ def getnlev(self): 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) +icon_nest = MODELIFNO(1377, 657, 51, 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 02f4ebe..ac02843 100644 --- a/src/plotlib.py +++ b/src/plotlib.py @@ -1,31 +1,21 @@ #!/usr/bin/python3 - -#this modul include all function to plot ICON Nest charts for eu and ce +from copy import deepcopy import numpy as np -#matplotlib +# matplotlib import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt -from matplotlib.colors import ListedColormap, BoundaryNorm, LinearSegmentedColormap -from matplotlib.cm import get_cmap -from matplotlib.ticker import MaxNLocator - -from mpl_toolkits.axes_grid1 import make_axes_locatable +from matplotlib.colors import LinearSegmentedColormap # ListedColormap, BoundaryNorm, -#cartopy +# cartopy import cartopy.crs as crs import cartopy.feature as cfeature from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER states_provinces = cfeature.NaturalEarthFeature( category='cultural', name='admin_0_boundary_lines_land', scale='10m', facecolor='none') -import os -import sys -import datetime -from copy import deepcopy - # own moduls import utilitylib as ut import kinematiclib as kin @@ -35,6 +25,7 @@ # ===================================================================================================================== # 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"] @@ -67,6 +58,7 @@ def eu_states(hour, start, projection=crs.EuroPP()): 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): 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) @@ -78,9 +70,9 @@ def ce_states(hour, start, projection=crs.EuroPP(), lon1=1.56, lon2=18.5, lat1=4 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): - # 10.77, 18.92, 49.80, 55.06 + + +def customize_area(hour, start, 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]) @@ -92,6 +84,7 @@ def east_de(hour, start, projection=crs.EuroPP(), lon1=10.7, lon2=18, lat1=49.8, 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): fig, ax = plt.subplots(figsize=(11, 9), subplot_kw=dict(projection=projection)) ax.set_extent([lon1, lon2, lat1, lat2]) @@ -117,10 +110,9 @@ def two_plots(projection=crs.EuroPP(), lon1=3.56, lon2=16.5, lat1=46.2, lat2=55. # --------------------------------------------------------------------------------------------------------------------- - # 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]) +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"]) # --------------------------------------------------------------------------------------------------------------------- @@ -144,7 +136,8 @@ def test_plot (cape_fld, lats, lons, hour, run, titel='CAPE'): fig, ax = ce_states(hour, run, projection=crs.PlateCarree()) 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) cax = fig.add_axes([0.27, 0.05, 0.35, 0.05]) fig.colorbar(wx, cax=cax, orientation='horizontal') @@ -193,19 +186,19 @@ 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 + # plot data wdir, spd = kin.uv2spddir(u, v) # smoothing - if smooth == True: + if smooth is True: wdir[1:] = (wdir[1:] + wdir[:-1])/2 spd[1:] = (spd[1:] + spd[:-1])/2 - # draw part of second cricle + # 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), - 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) @@ -213,7 +206,7 @@ def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): 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', threshold=10., imfmt="png"): +def basic_plot(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold=10., imfmt="png"): """ Parameters: ------------ @@ -233,19 +226,18 @@ def basic_plot (cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold 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) - + extend='max', alpha=0.4, antialiased=True) + 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]) > threshold: - hodopoint((lons[i, j], lats[i, j]), - np.mean(u[::-1, i-1:i+1, j-1:j+1], axis=(1, 2)), + 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') - - + 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) @@ -256,46 +248,48 @@ def basic_plot (cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold plt.savefig(name) plt.close() -def basic_plot_eastde (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"): - - fig, ax = east_de(hour, start, projection=crs.PlateCarree()) + 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) - + 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): 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()) - + 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') - - + 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-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"./images/hodographmap_east_{hour}.{imfmt}" + name = f"./images/hodographmap_area_{hour}.{imfmt}" plt.savefig(name) plt.close() # --------------------------------------------------------------------------------------------------------------------- -def nixon_hodograph(point, u, v, p, height, ax, width=0.1, clim=40, proj=crs.EuroPP()): +def nixon_hodograph(point, u, v, p, height, ax, width=0.1, clim=40, proj='polar', smooth=False): """ u, v : horizontal wind components rstu, rstv : storm motion vector for right mover """ - + i=0 - while height[i] < 500 : + while height[i] < 500: i += 1 i_500m = deepcopy(i) @@ -305,41 +299,37 @@ def nixon_hodograph(point, u, v, p, height, ax, width=0.1, clim=40, proj=crs.Eur while height[i] < 6000: i += 1 i_6km = deepcopy(i) - + rstu, rstv, lstu, lstv, mwu6, mwv6 = kin.non_parcel_bunkers_motion_experimental(u, v, p, i_500m, i_5km, i_6km) - + u-=rstu v-=rstv - - #plot + + # plot 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='polar') + 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.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) - #plot data + + # 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 is True: + wdir[1:] = (wdir[1:] + wdir[:-1])/2 + spd[1:] = (spd[1:] + spd[:-1])/2 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) @@ -350,7 +340,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, p, high, lats, lons, hour, start, imfmt="png"): +def nixon_proj(cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt="png"): """ Nixon projection cape_fld : 2D cape field @@ -370,28 +360,27 @@ def nixon_proj (cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt 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]) + # 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"]) - 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) - #cb.set_label(r'$m^2/s^2$') + 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) + # cb.set_label(r'$m^2/s^2$') - #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 + # 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): - #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)), + 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]) fig.colorbar(wx, cax=cax, orientation='horizontal') ax.annotate(r'$J/kg$', xy=(0.65, -0.04), xycoords='axes fraction', fontsize=14) diff --git a/src/utilitylib.py b/src/utilitylib.py index feeda9f..c743978 100644 --- a/src/utilitylib.py +++ b/src/utilitylib.py @@ -11,13 +11,12 @@ def datum (h, start): - z = '36' - a = type(z) - if (type(h) == a): + if not isinstance(h, int): h = int(h) - # type(h) - if (type(start) == a): + + 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) @@ -45,6 +44,7 @@ def load_yaml(yaml_file, yaml_path='./'): # 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="./"): opendataserver="https://opendata.dwd.de/weather/nwp/icon-eu/grib" @@ -57,9 +57,9 @@ def download_nwp(fieldname, datum="20240227", run="00", fp=0, store_path="./"): print("Download complete.") - # --------------------------------------------------------------------------------------------------------------------- + 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" @@ -95,25 +95,24 @@ def open_gribfile_multi(fieldname, lvl, datetime_obj, run, fp, path="./iconnest/ # 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/"): - data = Dataset(path+fieldname, 'r') + data = Dataset(f"{path}{fieldname}", 'r') lat = data.variables['lat'][:] lon = data.variables['lon'][:] - lons, lats= np.meshgrid(lon, lat) + lons, lats = np.meshgrid(lon, lat) - data = data.variables[fieldname.lower()][0,:,:,:].filled(np.nan) + data = data.variables[fieldname.lower()][0, :, :, :].filled(np.nan) data.close()