diff --git a/.gitignore b/.gitignore index 91024c3..a022747 100644 --- a/.gitignore +++ b/.gitignore @@ -1,7 +1,9 @@ # ignore ALL .log files -.ipynb_checkpoints/ -*/.ipynb_checkpoints/ +**__pycache__ + +**.ipynb_checkpoints/ + */*.png */*/*.png diff --git a/README.md b/README.md index 8f1f2ce..c47b923 100644 --- a/README.md +++ b/README.md @@ -16,9 +16,9 @@ You need to addg Conda initialization to your etc/profile, as well. ![Example](images/example.png) -![Map of Hodographs](images/hodographmap_area.png) +![Map of Hodographs IFS](images/hodographmap_IFS.png) -![Map of Soundings](images/soundingmap_example.png) +![Map of Hodographs GFS](images/hodographmap_GFS.png) ## How to run @@ -35,16 +35,16 @@ python3 main.py ICON cd images ``` +## Datasource +- [ICON Nest](https://opendata.dwd.de/weather/nwp/icon-eu/) +- [IFS](https://www.ecmwf.int/en/forecasts/datasets/open-data) +- [GFS (NOAA)](https://www.nco.ncep.noaa.gov/pmb/products/gfs/) + ## 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 -- GFS: https://www.nco.ncep.noaa.gov/pmb/products/gfs/ - ## License This project is licensed under the terms of the MIT license. See the [LICENSE](LICENSE) file for details. diff --git a/images/hodographmap_GFS.png b/images/hodographmap_GFS.png new file mode 100644 index 0000000..b11732e Binary files /dev/null and b/images/hodographmap_GFS.png differ diff --git a/images/hodographmap_IFS.png b/images/hodographmap_IFS.png new file mode 100644 index 0000000..8050b25 Binary files /dev/null and b/images/hodographmap_IFS.png differ diff --git a/src/__pycache__/kinematiclib.cpython-311.pyc b/src/__pycache__/kinematiclib.cpython-311.pyc deleted file mode 100644 index 33e17b3..0000000 Binary files a/src/__pycache__/kinematiclib.cpython-311.pyc and /dev/null differ diff --git a/src/main.py b/src/main.py index f7c6180..a9b75f5 100644 --- a/src/main.py +++ b/src/main.py @@ -15,23 +15,29 @@ # --------------------------------------------------------------------------------------------------------------------- -def run(model_obj, program_mode, fieldname, rundate, model_run, fp): +def run(model_obj, fp): config = ut.load_yaml('config.yml') - print(model_obj) + if model_obj.getname() == "IFS" or model_obj.getname() == "GFS": - cape_fld, u_fld, v_fld, lats, lons = ut.open_gribfile_preslvl(model_obj, rundate, model_run, fp, path="./modeldata/") - pres_levels = model_obj.getlevels() + cape_fld, u_fld, v_fld, lats, lons = ut.open_gribfile_preslvl(model_obj, fp, path="./modeldata/") + print(np.nanmean(u_fld, axis=(1, 2))) - plotlib.basic_plot(cape_fld, u_fld, v_fld, pres_levels, 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"]) else: + # program_mode + program_mode = config["program_mode"] + + # replace space with underscores + fieldname = "CAPE_ML" #args.field.replace(" ", "_") + rundate = model_obj.getrundate() if program_mode == "Test": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./modeldata/") + cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_obj.getrun(), fp, path="./modeldata/") assert cape_fld.shape == (model_obj.getnlat(), model_obj.getnlon()), "Shape inconsistency" - plotlib.test_plot(cape_fld, lats, lons, fp, model_run, titel='CAPE') + plotlib.test_plot(cape_fld, lats, lons, fp, model_obj.getrun(), titel='CAPE') elif program_mode == "Basic": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./modeldata/") + cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_obj.getrun(), fp, path="./modeldata/") nlvl = int(model_obj.getnlev()) u_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) @@ -42,18 +48,18 @@ def run(model_obj, program_mode, fieldname, rundate, model_run, fp): lvl_idx = 0 pres_levels = model_obj.getlevels() for level in pres_levels: - u_fld[lvl_idx, :, :] = ut.open_icon_gribfile_preslvl("U", level, rundate, model_run, fp, path="./modeldata/") - v_fld[lvl_idx, :, :] = ut.open_icon_gribfile_preslvl("V", level, rundate, model_run, fp, path="./modeldata/") + u_fld[lvl_idx, :, :] = ut.open_icon_gribfile_preslvl("U", level, rundate, model_obj.getrun(), fp, path="./modeldata/") + v_fld[lvl_idx, :, :] = ut.open_icon_gribfile_preslvl("V", level, rundate, model_obj.getrun(), 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, pres_levels, lats, lons, fp, model_run, - titel='CAPE', threshold=config["threshold"]) - plotlib.basic_plot_custarea(cape_fld, u_fld, v_fld, pres_levels, 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"]) else: print("Wrong command line argument") exit(-1) @@ -62,7 +68,6 @@ def run(model_obj, program_mode, fieldname, rundate, model_run, fp): def main(): - config = ut.load_yaml('config.yml') run_config = ut.load_yaml('run.yml') # command line arguments @@ -94,21 +99,6 @@ def main(): print("Unknown field") exit(-1) - if args.date is None: - rundate = datetime.strptime(run_config["default_date"], "%Y-%m-%d") - else: - rundate = datetime.strptime(args.date, "%Y-%m-%d") - - if args.fp is None: - fp = run_config["fp"] - else: - fp = args.fp - - if args.run is None: - model_run = run_config["run"] - else: - model_run = args.run - if args.Model is None: model_obj = model.MODELIFNO("ICON EU", 1377, 657, 0.0625, "pres") elif "ICON" in args.Model: @@ -121,15 +111,24 @@ def main(): print("Unkown model! Exit.") exit(0) - # program_mode - program_mode = config["program_mode"] + if args.date is None: + model_obj.setrundate(datetime.strptime(run_config["default_date"], "%Y-%m-%d")) + else: + model_obj.setrundate(datetime.strptime(args.date, "%Y-%m-%d")) - # replace space with underscores - fieldname = args.field.replace(" ", "_") + if args.fp is None: + fp = run_config["fp"] + else: + fp = args.fp + + if args.run is None: + model_obj.setrun(run_config["run"]) + else: + model_obj.setrun(args.run) - print(f"\nDate: {rundate}\n Arguments: {args} \nConfig-File: {config}\n\n") + print(f"\nArguments: {args}\n{model_obj}") - run(model_obj, program_mode, fieldname, rundate, model_run, fp) + run(model_obj, fp) # --------------------------------------------------------------------------------------------------------------------- diff --git a/src/meteolib.py b/src/meteolib.py index 83f8778..10187e7 100644 --- a/src/meteolib.py +++ b/src/meteolib.py @@ -21,10 +21,9 @@ def uv2spddir(u, v): direction = np.rad2deg(np.arctan2(-u, v)) if isinstance(direction, np.ndarray): - direction[np.where(direction < 0)] += 360 + direction = np.remainder(direction + 360, 360) else: - if direction < 0: - direction += 360 + direction = (direction + 360) % 360 return (np.deg2rad(direction), np.sqrt(u*u + v*v)) diff --git a/src/modelinfolib.py b/src/modelinfolib.py index 1c35a8c..2716e8c 100644 --- a/src/modelinfolib.py +++ b/src/modelinfolib.py @@ -1,4 +1,7 @@ #!/usr/bin/python3 +from datetime import datetime, date +from utilitylib import load_yaml + """ # ICON Nest points = 904689 @@ -51,8 +54,6 @@ ('u', 'isobaricInhPa', 500), ('u', 'isobaricInhPa', 400), ('u', 'isobaricInhPa', 300), - ('u', 'isobaricInhPa', 250), - ('u', 'isobaricInhPa', 200), ('v', 'isobaricInhPa', 1000), ('v', 'isobaricInhPa', 975), ('v', 'isobaricInhPa', 950), @@ -67,8 +68,6 @@ ('v', 'isobaricInhPa', 500), ('v', 'isobaricInhPa', 400), ('v', 'isobaricInhPa', 300), - ('v', 'isobaricInhPa', 250), - ('v', 'isobaricInhPa', 200), ('cape', 'pressureFromGroundLayer', 9000), # 9000 18000 25500 # ('cape', 'surface', 0), # ('cape', 'pressureFromGroundLayer', 9000), # 9000 18000 25500 @@ -86,8 +85,6 @@ ('u', 'isobaricInhPa', 500), ('u', 'isobaricInhPa', 400), ('u', 'isobaricInhPa', 300), - ('u', 'isobaricInhPa', 250), - ('u', 'isobaricInhPa', 200), ('v', 'isobaricInhPa', 1000), ('v', 'isobaricInhPa', 925), ('v', 'isobaricInhPa', 850), @@ -96,8 +93,6 @@ ('v', 'isobaricInhPa', 500), ('v', 'isobaricInhPa', 400), ('v', 'isobaricInhPa', 300), - ('v', 'isobaricInhPa', 250), - ('v', 'isobaricInhPa', 200), ('cape', 'entireAtmosphere', 0), # ('10u', 'heightAboveGround', 10), # ('10v', 'heightAboveGround', 10), @@ -111,32 +106,63 @@ class MODELIFNO: def __init__(self, modelname, nlon, nlat, d_grad, levtyp): + config = load_yaml('config.yml') self.modelname = modelname self.points = nlon*nlat self.nlon = nlon self.nlat = nlat if "ICON" in modelname: - self.levels = [1000, 950, 925, 900, 875, 850, 825, 800, 775, 700, 600, 500, 400, 300, 250, 200] + self.levels = [1000, 950, 925, 900, 875, 850, 825, 800, 775, 700, 600, 500, 400, 300] self.parlist = None + self.hodo_interval_lat = range(272, 415, 12) + self.hodo_interval_lon = range(420, 670, 15) elif modelname == "IFS": - self.levels = [1000, 925, 850, 700, 600, 500, 400, 300, 250, 200] + self.levels = [1000, 925, 850, 700, 600, 500, 400, 300] self.parlist = par_list_ifs + self.hodo_interval_lat = range(143, 174, 3) + self.hodo_interval_lon = range(731, 794, 3) elif modelname == "GFS": - self.levels = [1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 500, 400, 300, 250, 200] + self.levels = [1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 500, 400, 300] self.parlist = par_list_gfs + self.hodo_interval_lat = range(143, 174, 3) + self.hodo_interval_lon = range(11, 74, 3) else: - self.levels = [1000, 850, 700, 600, 500, 400, 300, 250, 200] + self.levels = [1000, 850, 700, 600, 500, 400, 300] self.d_grad = d_grad self.levtyp = levtyp + self.run = -99 + self.rundate = date.today() def __str__(self): return ( - f"Model Information: {self.modelname}\nPoints: {self.points}\n" - f"Number of Longitudes: {self.nlon}\nNumber of Latitudes: {self.nlat}\n" + f"Model Information: {self.modelname} Run: {self.run} Date: {self.rundate}\n" + f"Points: {self.points}\n" + f"Number of Longitudes: {self.nlon}\tNumber of Latitudes: {self.nlat}\n" f"Horizontal resolution: {self.d_grad}\n" f"Number of Levels: {len(self.levels)}\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 @@ -170,8 +196,10 @@ 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("ICON EU", 1377, 657, 0.0625, "pres") +icon_nest = MODELIFNO("ICON Nest", 1377, 657, 0.0625, "pres") ifs = MODELIFNO("IFS", 1440, 721, 0.25, "pres") -gfs = MODELIFNO("GFS", 1440, 721, 0.25, "pres") \ No newline at end of file +gfs = MODELIFNO("GFS", 1440, 721, 0.25, "pres") diff --git a/src/plotlib.py b/src/plotlib.py index 738b34b..18dedee 100644 --- a/src/plotlib.py +++ b/src/plotlib.py @@ -20,7 +20,6 @@ # own moduls import utilitylib as ut import meteolib as met -from skewT import SkewXAxes # --------------------------------------------------------------------------------------------------------------------- # create plot class @@ -33,7 +32,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 +42,60 @@ 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, model_name, 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(model_name, xy=(0.02, -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 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 @@ -190,7 +189,6 @@ def hodopoint(point, u, v, pres_levels, ax, width=0.1, clim=40, proj='polar', sm 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 @@ -199,10 +197,14 @@ def hodopoint(point, u, v, pres_levels, ax, width=0.1, clim=40, proj='polar', sm spd[1:] = (spd[1:] + spd[:-1])/2 # draw part of second cricle - if np.max(spd) > 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) + if np.max(spd[4:]) > 28: + u_mean = np.mean(u[np.where(spd > 27)]) + v_mean = np.mean(v[np.where(spd > 27)]) + 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) + # plot data idx_low = pres_levels.index(850) idx_mid = pres_levels.index(600) ax2.plot(wdir[:idx_low+1:1], spd[:idx_low+1:1], 'r-', lw=1.5) @@ -213,7 +215,7 @@ def hodopoint(point, u, v, pres_levels, ax, width=0.1, clim=40, proj='polar', sm # --------------------------------------------------------------------------------------------------------------------- -def basic_plot(cape_fld, u, v, pres_levels, 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: ------------ @@ -229,45 +231,50 @@ def basic_plot(cape_fld, u, v, pres_levels, lats, lons, hour, start, titel='CAPE None """ - fig, ax = ce_states(hour, start, projection=crs.PlateCarree()) - plt.title(titel, fontsize=titlesize) + pres_levels = model_obj.getlevels() + 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=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 range(142, 176, 4): - for j in range(731, 794, 4): - if np.mean(cape_fld[i-1:i+1, j-1:j+1]) > 0: + 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)), pres_levels, 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: srf-850hPa', xy=(0.75, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('green: 850-600hPa', xy=(0.75, -0.07), xycoords='axes fraction', fontsize=14) - ax.annotate('blue: above 600hPa', 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('1000-850 hPa in red', xy=(0.02, -0.03), xycoords='axes fraction', fontsize=13) + ax.annotate(' 850-600 hPa in green', xy=(0.02, -0.06), xycoords='axes fraction', fontsize=13) + ax.annotate(' 600-250 hPa in blue', xy=(0.02, -0.09), xycoords='axes fraction', fontsize=13) + ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.12), xycoords='axes fraction', fontsize=13) - 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, pres_levels, 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"): + pres_levels = model_obj.getlevels() + model_name = model_obj.getname() + start = model_obj.getrun() 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) + fig, ax = customize_area(hour, start, model_obj.getrundate(), model_name, + 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) @@ -282,129 +289,13 @@ def basic_plot_custarea(cape_fld, u, v, pres_levels, lats, lons, hour, start, ti cax = fig.add_axes([0.27, 0.05, 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(r'$J/kg$', xy=(0.65, -0.04), xycoords='axes fraction', fontsize=14) ax.annotate('red: srf-850hPa', xy=(0.75, -0.04), xycoords='axes fraction', fontsize=14) ax.annotate('green: 850-600hPa', xy=(0.75, -0.07), xycoords='axes fraction', fontsize=14) ax.annotate('blue: above 600hPa', 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_area_{hour}.{imfmt}" - plt.savefig(name) - plt.close() - -# --------------------------------------------------------------------------------------------------------------------- - - -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: - i += 1 - i_500m = deepcopy(i) - - while height[i] < 5500: - i += 1 - i_5km = deepcopy(i) - while height[i] < 6000: - i += 1 - i_6km = deepcopy(i) - - 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 - - # plot - test = ax.transData.transform(point) - # this should take us from the display coordinates to the axes coordinates. - trans = ax.transAxes.inverted().transform(test) - 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_ylim(0, clim) - ax2.set_theta_offset(np.pi/2) - # 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) - - # plot data - wdir, spd = met.uv2spddir(u, v) - - # 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) - 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) - - -def nixon_proj(cape_fld, dls_fld, u, v, p, high, lats, lons, hour, start, imfmt="png"): - """ - Nixon projection - cape_fld : 2D cape field - dls_fld : 2D deep layer shear field or brn shear ... - u, v : wind components - high : model level high - or - rstu, rstv : storm motion vector - background filed is cape ... - only hodographs with more the 10 m/s dls - """ - - titel = 'CAPE with Hodographs' - - fig, ax = ce_states(hour, start, projection=crs.PlateCarree()) - 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 = 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$') - - # 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 - 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) - - 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/nixon_ce_{hour}.{imfmt}" + name = f"./images/hodographmap_area_{model_name.replace(' ', '_')}_{hour}.{imfmt}" plt.savefig(name) plt.close() diff --git a/src/skewT.py b/src/skewT.py deleted file mode 100644 index eebb494..0000000 --- a/src/skewT.py +++ /dev/null @@ -1,162 +0,0 @@ -from matplotlib.axes import Axes -import matplotlib.transforms as transforms -import matplotlib.axis as maxis -import matplotlib.spines as mspines -# from matplotlib.projections import register_projection - - -# The sole purpose of this class is to look at the upper, lower, or total -# interval as appropriate and see what parts of the tick to draw, if any. -class SkewXTick(maxis.XTick): - def update_position(self, loc): - # This ensures that the new value of the location is set before - # any other updates take place - self._loc = loc - super().update_position(loc) - - def _has_default_loc(self): - return self.get_loc() is None - - def _need_lower(self): - return (self._has_default_loc() or - transforms.interval_contains(self.axes.lower_xlim, - self.get_loc())) - - def _need_upper(self): - return (self._has_default_loc() or - transforms.interval_contains(self.axes.upper_xlim, - self.get_loc())) - - @property - def gridOn(self): - return (self._gridOn and (self._has_default_loc() or - transforms.interval_contains(self.get_view_interval(), - self.get_loc()))) - - @gridOn.setter - def gridOn(self, value): - self._gridOn = value - - @property - def tick1On(self): - return self._tick1On and self._need_lower() - - @tick1On.setter - def tick1On(self, value): - self._tick1On = value - - @property - def label1On(self): - return self._label1On and self._need_lower() - - @label1On.setter - def label1On(self, value): - self._label1On = value - - @property - def tick2On(self): - return self._tick2On and self._need_upper() - - @tick2On.setter - def tick2On(self, value): - self._tick2On = value - - @property - def label2On(self): - return self._label2On and self._need_upper() - - @label2On.setter - def label2On(self, value): - self._label2On = value - - def get_view_interval(self): - return self.axes.xaxis.get_view_interval() - - -# This class exists to provide two separate sets of intervals to the tick, -# as well as create instances of the custom tick -class SkewXAxis(maxis.XAxis): - def _get_tick(self, major): - return SkewXTick(self.axes, None, '', major=major) - - def get_view_interval(self): - return self.axes.upper_xlim[0], self.axes.lower_xlim[1] - - -# This class exists to calculate the separate data range of the -# upper X-axis and draw the spine there. It also provides this range -# to the X-axis artist for ticking and gridlines -class SkewSpine(mspines.Spine): - def _adjust_location(self): - pts = self._path.vertices - if self.spine_type == 'top': - pts[:, 0] = self.axes.upper_xlim - else: - pts[:, 0] = self.axes.lower_xlim - - -# This class handles registration of the skew-xaxes as a projection as well -# as setting up the appropriate transformations. It also overrides standard -# spines and axes instances as appropriate. -class SkewXAxes(Axes): - # The projection must specify a name. This will be used be the - # user to select the projection, i.e. ``subplot(111, - # projection='skewx')``. - name = 'skewx' - - def _init_axis(self): - # Taken from Axes and modified to use our modified X-axis - self.xaxis = SkewXAxis(self) - self.spines['top'].register_axis(self.xaxis) - self.spines['bottom'].register_axis(self.xaxis) - self.yaxis = maxis.YAxis(self) - self.spines['left'].register_axis(self.yaxis) - self.spines['right'].register_axis(self.yaxis) - - def _gen_axes_spines(self): - spines = {'top': SkewSpine.linear_spine(self, 'top'), - 'bottom': mspines.Spine.linear_spine(self, 'bottom'), - 'left': mspines.Spine.linear_spine(self, 'left'), - 'right': mspines.Spine.linear_spine(self, 'right')} - return spines - - def _set_lim_and_transforms(self): - """ - This is called once when the plot is created to set up all the - transforms for the data, text and grids. - """ - rot = 45 - - # Get the standard transform setup from the Axes base class - Axes._set_lim_and_transforms(self) - - # Need to put the skew in the middle, after the scale and limits, - # but before the transAxes. This way, the skew is done in Axes - # coordinates thus performing the transform around the proper origin - # We keep the pre-transAxes transform around for other users, like the - # spines for finding bounds - self.transDataToAxes = self.transScale + \ - self.transLimits + transforms.Affine2D().skew_deg(rot, 0) - - # Create the full transform from Data to Pixels - self.transData = self.transDataToAxes + self.transAxes - - # Blended transforms like this need to have the skewing applied using - # both axes, in axes coords like before. - self._xaxis_transform = (transforms.blended_transform_factory( - self.transScale + self.transLimits, - transforms.IdentityTransform()) + - transforms.Affine2D().skew_deg(rot, 0)) + self.transAxes - - @property - def lower_xlim(self): - return self.axes.viewLim.intervalx - - @property - def upper_xlim(self): - pts = [[0., 1.], [1., 1.]] - return self.transDataToAxes.inverted().transform(pts)[:, 0] - - -# Now register the projection with matplotlib so the user can select -# it. diff --git a/src/utilitylib.py b/src/utilitylib.py index 62471d3..a33c216 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 # --------------------------------------------------------------------------------------------------------------------- @@ -100,8 +98,8 @@ def open_icon_gribfile_preslvl(fieldname, lvl, datetime_obj, run, fp, path="./mo print("Maximum:", np.nanmax(data)) return data -def open_gribfile_preslvl(model_obj, datetime_obj, run, fp, path="./modeldata/"): - date_string = datetime_obj.strftime("%Y%m%d") +def open_gribfile_preslvl(model_obj, fp, path="./modeldata/"): + date_string = model_obj.getrundate_as_str("%Y%m%d") # create numpy.array fields shape = (model_obj.getnlev(), model_obj.getnlat(), model_obj.getnlon()) @@ -113,17 +111,17 @@ def open_gribfile_preslvl(model_obj, datetime_obj, run, fp, path="./modeldata/") lats = np.full(shape, np.nan) lons = np.full(shape, np.nan) pres_levels = model_obj.getlevels() + run = model_obj.getrun() # Open the GRIB file # modelname_RRz_YYYYMMDD_f015.grib2 - gribidx = pygrib.index(f"{path}{model_obj.getname().lower()}_{run:02d}z_{date_string}_f{fp:03d}.grib2",'shortName','typeOfLevel','level') - #grbs.seek(0) - + gribidx = pygrib.index(f"{path}{model_obj.getname().lower()}_{run:02d}z_{date_string}_f{fp:03d}.grib2", + 'shortName', 'typeOfLevel', 'level') + # grbs.seek(0) for par in model_obj.getParamter(): try: - grb_message = gribidx.select(shortName=par[0],typeOfLevel=par[1],level=par[2])[0] # takes the matching grib message - print(grb_message) + grb_message = gribidx.select(shortName=par[0], typeOfLevel=par[1], level=par[2])[0] # takes the matching grib message if par[0] == "cape": cape_fld = grb_message.values lats, lons = grb_message.latlons()