diff --git a/.github/workflows/python-nose.yml b/.github/workflows/python-nose.yml new file mode 100644 index 0000000..06d295b --- /dev/null +++ b/.github/workflows/python-nose.yml @@ -0,0 +1,24 @@ +name: Unittest + +on: + push: + branches: [ "main", "modellevel", "preslevel", "lint" ] + +jobs: + build-linux: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v3 + - name: Set up Python 3.11 + uses: actions/setup-python@v3 + with: + python-version: '3.11' + - name: Install dependencies + run: | + pip3 install --upgrade pip + pip3 install -r requirements.txt + - name: Run tests with nose2 + run: | + cd test + nose2 diff --git a/.gitignore b/.gitignore index 96cb0ae..a022747 100644 --- a/.gitignore +++ b/.gitignore @@ -1,7 +1,9 @@ # ignore ALL .log files -.ipynb_checkpoints/ -*/.ipynb_checkpoints/ +**__pycache__ + +**.ipynb_checkpoints/ + */*.png */*/*.png @@ -16,5 +18,4 @@ */*.pyc */*/*.pyc -src/invariant -src/iconnest +src/modeldata diff --git a/README.md b/README.md index e830020..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 @@ -29,24 +29,22 @@ bash download_script.sh conda activate HodographMaps # Plot Hodograph -python3 main.py Basic - -# Plot Soundings -python3 main.py Sounding - +python3 main.py IFS +python3 main.py GFS +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/environment.yml b/environment.yml index c13e6bc..592a723 100644 --- a/environment.yml +++ b/environment.yml @@ -10,3 +10,4 @@ dependencies: - pygrib - netcdf4 - pyyaml + - nose2=0.9.2 diff --git a/images/example.png b/images/example.png index 8f42aa7..9d88cc5 100644 Binary files a/images/example.png and b/images/example.png differ 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/images/hodographmap_area.png b/images/hodographmap_area.png index 2572ad0..6d1b4e7 100644 Binary files a/images/hodographmap_area.png and b/images/hodographmap_area.png differ diff --git a/images/soundingmap_example.png b/images/soundingmap_example.png deleted file mode 100644 index 60043dd..0000000 Binary files a/images/soundingmap_example.png and /dev/null differ diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..4e72581 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,8 @@ +numpy +matplotlib +cartopy +requests +pygrib +netcdf4 +pyyaml +nose2==0.9.2 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/config.yml b/src/config.yml index 0b03f7b..9f78272 100644 --- a/src/config.yml +++ b/src/config.yml @@ -1,8 +1,6 @@ -levels: [74,20] -steps: 1 # level steps +program_mode : Basic # Test Basic hodomap : True -soundmap : True threshold: 0 @@ -16,4 +14,3 @@ customize: # east de lat1: 50.1 lat2: 54.8 -# 10.77, 18.92, 49.80, 55.06 diff --git a/src/download_script.sh b/src/download_script.sh index 116fc04..e3746db 100755 --- a/src/download_script.sh +++ b/src/download_script.sh @@ -4,7 +4,7 @@ source /etc/profile conda activate HodographMaps -store_path=$(pwd)/iconnest +store_path=$(pwd)/modeldata # create nwp directory and if not there a output images directory mkdir -p ${store_path} @@ -12,24 +12,33 @@ mkdir -p ./images # select run R=0 +R=15 + # Path of the icon nest on the opendata-sever -model_pfad=https://opendata.dwd.de/weather/nwp/icon-eu/grib/$(printf "%02d" "$R") +icon_model_pfad=https://opendata.dwd.de/weather/nwp/icon-eu/grib/$(printf "%02d" "$R") +# Path the ifs on the opendata-sever +ifs_model_pfad=https://data.ecmwf.int/forecasts #/20240406/00z/ifs/0p4-beta/oper/ +# Path the gfs on the opendata-sever +gfs_model_pfad=https://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod + # date D=$(date +"%Y%m%d") +D=20240415 -echo "ICON EU NEST Run: " ${R} +echo "Run: " ${R} echo "Date: " $(date) echo "--------------------------------" ####################################################################### # name of the model files -regular_single=icon-eu_europe_regular-lat-lon_single-level_${D}$(printf "%02d" "$R")_ -regular_pressure=icon-eu_europe_regular-lat-lon_pressure-level_${D}$(printf "%02d" "$R")_ -regular_model=icon-eu_europe_regular-lat-lon_model-level_${D}$(printf "%02d" "$R")_ -invariant=icon-eu_europe_regular-lat-lon_time-invariant_${D}$(printf "%02d" "$R")_ +icon_single=icon-eu_europe_regular-lat-lon_single-level_${D}$(printf "%02d" "$R")_ +icon_pressure=icon-eu_europe_regular-lat-lon_pressure-level_${D}$(printf "%02d" "$R")_ + + + -for X in 15 #9 12 15 18 21 24 27 30 33 36 39 42 45 48 51 54 57 60 +for X in 3 #9 12 15 18 21 24 27 30 33 36 39 42 45 48 51 54 57 60 do T=$(printf "%03d" "$X") # single level @@ -37,21 +46,34 @@ do do typeset -l nvar nvar=${N} - wget -q ${model_pfad}/${nvar}/${regular_single}${T}_${N}.grib2.bz2 -P ${store_path} - bzip2 -dfq ${store_path}/${regular_single}${T}_${N}.grib2.bz2 + wget -q ${icon_model_pfad}/${nvar}/${icon_single}${T}_${N}.grib2.bz2 -P ${store_path} + bzip2 -dfq ${store_path}/${icon_single}${T}_${N}.grib2.bz2 done - for H in {20..74} + for H in 1000 950 925 900 875 850 825 800 775 700 600 500 400 300 250 200 do - for N in U V T QV P + for N in U V do typeset -l nvar nvar=${N} - wget -q ${model_pfad}/${nvar}/${regular_model}${T}_${H}_${N}.grib2.bz2 -P ${store_path} - bzip2 -dfq ${store_path}/${regular_model}${T}_${H}_${N}.grib2.bz2 + wget -q ${icon_model_pfad}/${nvar}/${icon_pressure}${T}_${H}_${N}.grib2.bz2 -P ${store_path} + bzip2 -dfq ${store_path}/${icon_pressure}${T}_${H}_${N}.grib2.bz2 done done + # ifs + ifs_file=${ifs_model_pfad}/${D}/$(printf "%02d" "$R")z/ifs/0p25/oper/${D}$(printf "%02d" "$R")0000-${X}h-oper-fc.grib2 + ifs_index=${ifs_model_pfad}/${D}/$(printf "%02d" "$R")z/ifs/0p25/oper/${D}$(printf "%02d" "$R")0000-${X}h-oper-fc.grib2 + wget ${ifs_file} -P ${store_path}/ + wget ${ifs_index} -P ${store_path}/ + mv ${store_path}/${D}$(printf "%02d" "$R")0000-${X}h-oper-fc.grib2 ${store_path}/ifs_$(printf "%02d" "$R")z_${D}_f${T}.grib2 + mv ${store_path}/${D}$(printf "%02d" "$R")0000-${X}h-oper-fc.index ${store_path}/ifs_$(printf "%02d" "$R")z_${D}_f${T}.index + + # gfs + gfs_file=${gfs_model_pfad}/gfs.${D}/$(printf "%02d" "$R")/atmos/gfs.t$(printf "%02d" "$R")z.pgrb2.0p25.f${T} + wget ${gfs_file} -P ${store_path}/ + mv ${store_path}/gfs.t$(printf "%02d" "$R")z.goessimpgrb2.0p25.f${T} ${store_path}/gfs_$(printf "%02d" "$R")z_${D}_f${T}.grib2 + # write run.yml echo run: ${R} > run.yml echo fp: ${X} >> run.yml @@ -65,4 +87,4 @@ do done # remove nwp files -rm -rf ${store_path} +#rm -rf ${store_path} diff --git a/src/main.py b/src/main.py index 292504a..35242db 100644 --- a/src/main.py +++ b/src/main.py @@ -11,132 +11,72 @@ import utilitylib as ut import plotlib import modelinfolib as model -model = model.icon_nest -# --------------------------------------------------------------------------------------------------------------------- +# --------------------------------------------------------------------------------------------------------------------------- -def run(program_mode, fieldname, rundate, model_run, fp): +def run(model_obj, fp): config = ut.load_yaml('config.yml') - if program_mode == "Test": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./iconnest/") - assert cape_fld.shape == (model.getnlat, model.getnlon), "Shape inconsistency" - plotlib.test_plot(cape_fld, lats, lons, fp, model_run, titel='CAPE') - - elif program_mode == "Sounding": - cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./iconnest/") - - 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) - q_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon())) - q_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): - 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, 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"] - 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, 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 + if model_obj.getname() == "IFS" or model_obj.getname() == "GFS": + 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, 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) - 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, 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))) - - du = np.subtract(u_fld[30, :, :], u_fld[0, :, :]) - dv = np.subtract(v_fld[30, :, :], v_fld[0, :, :]) - dls_fld = np.sqrt(np.add(np.square(du), np.square(dv))) - plotlib.nixon_proj(cape_fld, dls_fld, u_fld, v_fld, p_fld, h_fld, lats, lons, fp, model_run, imfmt="png") + plotlib.basic_plot(model_obj, cape_fld, u_fld, v_fld, lats, lons, fp, + threshold=config["threshold"]) else: - print("Wrong command line argument") - exit(-1) - -# --------------------------------------------------------------------------------------------------------------------- + # 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_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_obj.getrun(), titel='CAPE') + + elif program_mode == "Basic": + 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())) + u_fld.fill(np.nan) + v_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon())) + v_fld.fill(np.nan) + + 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_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(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) + +# --------------------------------------------------------------------------------------------------------------------------- def main(): - config = ut.load_yaml('config.yml') run_config = ut.load_yaml('run.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') + parser.add_argument('Model', type=str, + help='Model: ICON, IFS, or GFS') # Add optional argument parser.add_argument('-f', '--field', type=str, @@ -161,10 +101,22 @@ def main(): print("Unknown field") exit(-1) + if args.Model is None: + model_obj = model.MODELIFNO("ICON EU", 1377, 657, 0.0625, "pres") + elif "ICON" in args.Model: + model_obj = model.icon_nest + elif args.Model == "IFS": + model_obj = model.ifs + elif args.Model == "GFS": + model_obj = model.gfs + else: + print("Unkown model! Exit.") + exit(0) + if args.date is None: - rundate = datetime.strptime(run_config["default_date"], "%Y-%m-%d") + model_obj.setrundate(datetime.strptime(run_config["default_date"], "%Y-%m-%d")) else: - rundate = datetime.strptime(args.date, "%Y-%m-%d") + model_obj.setrundate(datetime.strptime(args.date, "%Y-%m-%d")) if args.fp is None: fp = run_config["fp"] @@ -172,22 +124,13 @@ def main(): fp = args.fp if args.run is None: - model_run = run_config["run"] - else: - model_run = args.run - - if args.mode != "Test" and args.mode != "Sounding" and args.mode != "Basic" and args.mode != "Nixon": - print("Unknown Mode. Exit program.") - exit(0) + model_obj.setrun(run_config["run"]) else: - program_mode = args.mode - - # replace space with underscores - fieldname = args.field.replace(" ", "_") + 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(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 67e3402..5b9043b 100644 --- a/src/modelinfolib.py +++ b/src/modelinfolib.py @@ -1,10 +1,12 @@ #!/usr/bin/python3 +from datetime import datetime, date + """ # ICON Nest points = 904689 nlon = 1377 nlat = 657 -nlev = 51 +lev = [1000, 950, 925, 900, 875, 850, 825, 800, 775, 700, 600, 500, 400, 300, 250, 200] lonmin = -23.5 lonmax = 45.0 latmin = 29.5 @@ -12,21 +14,21 @@ d_grad = 0.0625 # IFS -points = 405900 -nlon = 900 -nlat = 451 -nlev = +points = 1038240 +nlon = 1440 +nlat = 721 +lev = [1000, 925, 850, 700, 600, 500, 400, 300, 250, 200] lonmin = 0 lonmax = 359 latmin = -90 latmax = 90 -d_grad = 0.4 +d_grad = 0.25 # GFS points = 1038240 nlon = 1440 nlat = 721 -nlev = +nlev = [1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 500, 400, 300, 250, 200] lonmin = 0 lonmax = 359 latmin = -90 @@ -34,36 +36,138 @@ d_grad = 0.25 """ +# --------------------------------------------------------------------------------------------------------------------- +# name typeOfLevel level +par_list_gfs = [ + ('u', 'isobaricInhPa', 1000), + ('u', 'isobaricInhPa', 975), + ('u', 'isobaricInhPa', 950), + ('u', 'isobaricInhPa', 925), + ('u', 'isobaricInhPa', 900), + ('u', 'isobaricInhPa', 850), + ('u', 'isobaricInhPa', 800), + ('u', 'isobaricInhPa', 750), + ('u', 'isobaricInhPa', 700), + ('u', 'isobaricInhPa', 650), + ('u', 'isobaricInhPa', 600), + ('u', 'isobaricInhPa', 500), + ('u', 'isobaricInhPa', 400), + ('u', 'isobaricInhPa', 300), + ('v', 'isobaricInhPa', 1000), + ('v', 'isobaricInhPa', 975), + ('v', 'isobaricInhPa', 950), + ('v', 'isobaricInhPa', 925), + ('v', 'isobaricInhPa', 900), + ('v', 'isobaricInhPa', 850), + ('v', 'isobaricInhPa', 800), + ('v', 'isobaricInhPa', 750), + ('v', 'isobaricInhPa', 700), + ('v', 'isobaricInhPa', 650), + ('v', 'isobaricInhPa', 600), + ('v', 'isobaricInhPa', 500), + ('v', 'isobaricInhPa', 400), + ('v', 'isobaricInhPa', 300), + ('cape', 'pressureFromGroundLayer', 9000), # 9000 18000 25500 + # ('cape', 'surface', 0), + # ('cape', 'pressureFromGroundLayer', 9000), # 9000 18000 25500 + # ('10u', 'heightAboveGround', 10), + # ('10v', 'heightAboveGround', 10), + # ('100u', 'heightAboveGround', 100), + # ('100v', 'heightAboveGround', 100), +] +par_list_ifs = [ + ('u', 'isobaricInhPa', 1000), + ('u', 'isobaricInhPa', 925), + ('u', 'isobaricInhPa', 850), + ('u', 'isobaricInhPa', 700), + ('u', 'isobaricInhPa', 600), + ('u', 'isobaricInhPa', 500), + ('u', 'isobaricInhPa', 400), + ('u', 'isobaricInhPa', 300), + ('v', 'isobaricInhPa', 1000), + ('v', 'isobaricInhPa', 925), + ('v', 'isobaricInhPa', 850), + ('v', 'isobaricInhPa', 700), + ('v', 'isobaricInhPa', 600), + ('v', 'isobaricInhPa', 500), + ('v', 'isobaricInhPa', 400), + ('v', 'isobaricInhPa', 300), + ('cape', 'entireAtmosphere', 0), + # ('10u', 'heightAboveGround', 10), + # ('10v', 'heightAboveGround', 10), + # ('100u', 'heightAboveGround', 100), + # ('100v', 'heightAboveGround', 100), +] + +# --------------------------------------------------------------------------------------------------------------------- class MODELIFNO: - def __init__(self, nlon, nlat, d_grad, nlev, levtyp): + 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 - self.nlev = nlev + if "ICON" in modelname: + 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] + 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] + 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] self.d_grad = d_grad self.levtyp = levtyp - - """ - def __init__(self, nlon, nlat, nlev, lonmin, lonmax, latmin, latmax, d_grad): - self.points = nlon*nlat - self.nlon = nlon - self.nlat = nlat - self.nlev = nlev - self.d_grad = d_grad - print(lonmin, lonmax, latmin, latmax) - """ + self.run = -99 + self.rundate = date.today() def __str__(self): return ( - f"Model Information:\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: {self.nlev}\tLeveltyps: {self.levtyp}\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 + + def getParamter(self): + return self.parlist + def getpoints(self): return self.points @@ -73,8 +177,17 @@ def getnlon(self): def getnlat(self): return self.nlat + def getlevels(self): + return self.levels + + def getpreslevel_by_idx(self, idx): + if idx >= 0 and idx < len(self.levels): + return self.levels[idx] + else: + return -99 + def getnlev(self): - return self.nlev + return len(self.levels) def getlevtyp(self): return self.d_grad @@ -82,6 +195,11 @@ def getlevtyp(self): def getd_grad(self): return self.d_grad + def create_plottitle(self): + return f"Hodographmap of {self.modelname}" + # Example usage: -icon_nest = MODELIFNO(1377, 657, 0.0625, 54, "model") # lowest 74 and we download till 20 +icon_nest = MODELIFNO("ICON Nest", 1377, 657, 0.0625, "pres") +ifs = MODELIFNO("IFS", 1440, 721, 0.25, "pres") +gfs = MODELIFNO("GFS", 1440, 721, 0.25, "pres") diff --git a/src/plotlib.py b/src/plotlib.py index 60351a3..346deac 100644 --- a/src/plotlib.py +++ b/src/plotlib.py @@ -1,26 +1,22 @@ #!/usr/bin/python3 -from copy import deepcopy - import numpy as np # matplotlib import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt -from matplotlib.colors import LinearSegmentedColormap # ListedColormap, BoundaryNorm, -from matplotlib.projections import register_projection +from matplotlib.colors import LinearSegmentedColormap, ListedColormap # BoundaryNorm, # 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') +states_provinces = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_boundary_lines_land', + scale='10m', facecolor='none') # own moduls import utilitylib as ut import meteolib as met -from skewT import SkewXAxes # --------------------------------------------------------------------------------------------------------------------- # create plot class @@ -33,7 +29,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 +39,59 @@ def eu_merc(hour, start, projection=crs.Mercator(), factor=3): gl.ylabels_right = False gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER - string1, string2 = ut.datum(hour, start) - plt.annotate(string1, xy=(0.8, 1), xycoords='axes fraction', fontsize=fontsize) - plt.annotate(string2, xy=(0, 1), xycoords='axes fraction', fontsize=fontsize) + string1, string2 = ut.datum(hour, start, datetime_obj) + plt.annotate(model_name, xy=(0.7, -0.04), xycoords='axes fraction', fontsize=10) + plt.annotate(string1, xy=(0.82, 1.01), xycoords='axes fraction', fontsize=fontsize) + plt.annotate(string2, xy=(0, 1.01), xycoords='axes fraction', fontsize=fontsize) return fig, ax -def eu_states(hour, start, projection=crs.EuroPP()): +def eu_states(hour, start, datetime_obj, model_name, projection=crs.EuroPP()): fig, ax = plt.subplots(figsize=(9, 7), subplot_kw=dict(projection=projection)) plt.subplots_adjust(left=0.05, right=0.97, bottom=0.05, top=0.95) ax.set_extent([-9.5, 33.0, 38.5, 58.5]) ax.coastlines('50m', linewidth=1.2) ax.add_feature(states_provinces, edgecolor='black') - string1, string2 = ut.datum(hour, start) - plt.annotate("ICON Nest (DWD)", xy=(0.7, -0.04), xycoords='axes fraction', fontsize=10) - plt.annotate(string1, xy=(0.8, 1), xycoords='axes fraction', fontsize=fontsize) - plt.annotate(string2, xy=(0, 1), xycoords='axes fraction', fontsize=fontsize) + string1, string2 = ut.datum(hour, start, datetime_obj) + plt.annotate(model_name, xy=(0.7, -0.04), xycoords='axes fraction', fontsize=10) + plt.annotate(string1, xy=(0.82, 1.01), xycoords='axes fraction', fontsize=fontsize) + plt.annotate(string2, xy=(0, 1.01), xycoords='axes fraction', fontsize=fontsize) return fig, ax -def ce_states(hour, start, projection=crs.EuroPP(), lon1=1.56, lon2=18.5, lat1=45.1, lat2=56.6): +def ce_states(hour, start, datetime_obj, projection=crs.EuroPP(), lon1=1.56, lon2=18.5, lat1=45.1, lat2=56.6): fig, ax = plt.subplots(figsize=(11, 9), subplot_kw=dict(projection=projection)) plt.subplots_adjust(left=0.05, right=0.99, bottom=0.05, top=0.95) ax.set_extent([lon1, lon2, lat1, lat2]) ax.coastlines('10m', linewidth=1.2) ax.add_feature(states_provinces, edgecolor='black') - string1, string2 = ut.datum(hour, start) - plt.annotate("ICON Nest (DWD)", xy=(0.02, -0.04), xycoords='axes fraction', fontsize=10) - plt.annotate(string1, xy=(0.835, 1.02), xycoords='axes fraction', fontsize=fontsize) - plt.annotate(string2, xy=(0, 1.02), xycoords='axes fraction', fontsize=fontsize) + string1, string2 = ut.datum(hour, start, datetime_obj) + plt.annotate(string1, xy=(0.82, 1.01), xycoords='axes fraction', fontsize=fontsize) + plt.annotate(string2, xy=(0, 1.01), xycoords='axes fraction', fontsize=fontsize) return fig, ax -def customize_area(hour, start, projection=crs.EuroPP(), lon1=10.7, lon2=18, lat1=49.8, lat2=54.8): +def customize_area(hour, start, datetime_obj, 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(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 @@ -117,9 +112,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]) -cmap = LinearSegmentedColormap.from_list("", ["green", "yellow", "orange", "red", "darkred", "darkmagenta"]) -cmap2 = LinearSegmentedColormap.from_list("", ["gold", "orange", "darkorange", "red", "darkred", "darkmagenta"]) - +# cmap = LinearSegmentedColormap.from_list("", ["green", "yellow", "orange", "red", "darkred", "darkmagenta"]) +cmap = LinearSegmentedColormap.from_list("", ["lightgoldenrodyellow", "orange", "red"]) +cmap = ListedColormap(cmap(np.linspace(0, 1, 256))[8:255]) # --------------------------------------------------------------------------------------------------------------------- @@ -131,7 +126,7 @@ def test_plot(cape_fld, lats, lons, hour, run, titel='CAPE'): lat : lon : hour : - run : + run : Returns: -------- @@ -154,103 +149,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: - ------------ - point : tuple of data coordinates (10, 55) - temp : temperature - width : width of added axes for the soundings as fig coordinates (examples: 0.5 are half of the image) - proj : 'polar' - - Returns: - -------- - None - """ - register_projection(SkewXAxes) - - # this takes us from the data coordinates to the display coordinates. - test = ax.transData.transform(point) - - # this should take us from the display coordinates to the axes coordinates. - trans = ax.transAxes.inverted().transform(test) - - # create new axes - ax2 = plt.axes([trans[0]-width/2, trans[1]-width/2, width, width]) # , projection='skewx') - - ax2.get_xaxis().set_visible(False) - ax2.get_yaxis().set_visible(False) - ax2.set_frame_on(False) - - # ax2.set_xlim(-, ) - # 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)+30, '-k', alpha=.3, lw=0.8) - - # calculate dewpoint - pres /= 100 - 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() - -# --------------------------------------------------------------------------------------------------------------------- - - -def sounding_plot(cape_fld, temp, q_fld, p_fld, lats, lons, hour, start, titel='CAPE', imfmt="png"): - """ - Parameters: - ------------ - cape_fld : CAPE field - temp : temperature - lats : - lons : - hour : - start : - - Returns: - -------- - None - """ - - 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, - 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() - - 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 dot line: freezing level', xy=(0.75, -0.1), xycoords='axes fraction', fontsize=14) - ax.annotate('blue lines: temperature', xy=(0.75, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('green lines: dewpoints', xy=(0.75, -0.07), xycoords='axes fraction', fontsize=14) - name = f"./images/soundingmap_ce_{hour}.{imfmt}" - plt.savefig(name) - plt.close() - -# --------------------------------------------------------------------------------------------------------------------- - - -def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): +def hodopoint(point, u, v, pres_levels, ax, width=0.1, clim=50, proj='polar', smooth=False): """ Parameters: ------------ @@ -286,7 +185,6 @@ def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): ax2.plot(np.linspace(0, 2*np.pi, 100), np.zeros(100)+10, '-k', alpha=.3, lw=0.8) # ax2.plot(np.linspace(0, 2*np.pi, 100), np.zeros(100)+30, '-k', alpha=.3, lw=0.8) - # plot data wdir, spd = met.uv2spddir(u, v) # smoothing @@ -295,20 +193,25 @@ def hodopoint(point, u, v, ax, width=0.1, clim=40, proj='polar', smooth=False): spd[1:] = (spd[1:] + spd[:-1])/2 # 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), + 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) - - 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) + # 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], '-', color='darkmagenta', lw=1.5) + ax2.plot(wdir[idx_low:idx_mid+1:1], spd[idx_low:idx_mid+1:1], '-', color='navy', lw=1.5) + ax2.plot(wdir[idx_mid:], spd[idx_mid:], '-', color='darkgreen', 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', threshold=10., imfmt="png"): +def basic_plot(model_obj, cape_fld, u, v, lats, lons, hour, threshold=10., imfmt="png"): """ Parameters: ------------ @@ -324,43 +227,50 @@ def basic_plot(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', threshold= 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 model_obj.hodo_interval_lat: + for j in model_obj.hodo_interval_lon: if np.mean(cape_fld[i-1:i+1, j-1:j+1]) > threshold: hodopoint((lons[i, j], lats[i, j]), np.mean(u[:, i-1:i+1, j-1:j+1], axis=(1, 2)), - np.mean(v[:, i-1:i+1, j-1:j+1], axis=(1, 2)), ax, width=0.1) # , proj=crs.PlateCarree() + np.mean(v[:, i-1:i+1, j-1:j+1], axis=(1, 2)), 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: 1-10 model level', xy=(0.75, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('green: 10-20 model level', xy=(0.75, -0.07), xycoords='axes fraction', fontsize=14) - ax.annotate('blue: 20-50 model level', xy=(0.75, -0.1), xycoords='axes fraction', fontsize=14) - ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.07), xycoords='axes fraction', fontsize=10) + ax.annotate("CAPE ML(contour plot)", xy=(0.8, -0.07), xycoords='axes fraction', fontsize=14) + ax.annotate(r'in $J/kg$', xy=(0.8, -0.1), xycoords='axes fraction', fontsize=14) + + ax.annotate('1000-850 hPa in magenta', xy=(0.02, -0.03), xycoords='axes fraction', fontsize=13) + ax.annotate(' 850-600 hPa in blue', xy=(0.02, -0.06), xycoords='axes fraction', fontsize=13) + ax.annotate(' 600-300 hPa in green', 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, 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) @@ -370,134 +280,19 @@ def basic_plot_custarea(cape_fld, u, v, lats, lons, hour, start, titel='CAPE', t if np.mean(cape_fld[i-1:i+1, j-1:j+1]) > threshold: hodopoint((lons[i, j], lats[i, j]), np.mean(u[:, i-1:i+1, j-1:j+1], axis=(1, 2)), - np.mean(v[:, i-1:i+1, j-1:j+1], axis=(1, 2)), ax, width=0.1) # , proj=crs.PlateCarree() + np.mean(v[:, i-1:i+1, j-1:j+1], axis=(1, 2)), 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: 1-10 model level', xy=(0.75, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('green: 10-20 model level', xy=(0.75, -0.07), xycoords='axes fraction', fontsize=14) - ax.annotate('blue: 20-50 model level', xy=(0.75, -0.1), xycoords='axes fraction', fontsize=14) - ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.07), xycoords='axes fraction', fontsize=10) - - 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("CAPE ML(contour plot)", xy=(0.8, -0.07), xycoords='axes fraction', fontsize=14) + ax.annotate(r'in $J/kg$', xy=(0.8, -0.1), xycoords='axes fraction', fontsize=14) - ax.annotate('red: 1-10 model level', xy=(0.75, -0.04), xycoords='axes fraction', fontsize=14) - ax.annotate('green: 10-20 model level', xy=(0.75, -0.07), xycoords='axes fraction', fontsize=14) - ax.annotate('blue: 20-40 model level', xy=(0.75, -0.1), xycoords='axes fraction', fontsize=14) - ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.07), xycoords='axes fraction', fontsize=10) + ax.annotate('1000-850 hPa in magenta', xy=(0.02, -0.025), xycoords='axes fraction', fontsize=13) + ax.annotate(' 850-600 hPa in blue', xy=(0.02, -0.05), xycoords='axes fraction', fontsize=13) + ax.annotate(' 600-300 hPa in green', xy=(0.02, -0.075), xycoords='axes fraction', fontsize=13) + ax.annotate("grey circles are 10 and 30m/s", xy=(0.02, -0.11), xycoords='axes fraction', fontsize=13) - 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/run.yml b/src/run.yml index 019cabf..f3d56c4 100644 --- a/src/run.yml +++ b/src/run.yml @@ -1,3 +1,3 @@ -run: 0 -fp: 15 -default_date: "2024-04-09" +run: 15 +fp: 3 +default_date: "2024-04-15" 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 ee4dae0..ebbae4f 100644 --- a/src/utilitylib.py +++ b/src/utilitylib.py @@ -1,6 +1,6 @@ #!/usr/bin/python3 -import datetime +from datetime import datetime, timedelta import requests import yaml import pygrib @@ -9,22 +9,20 @@ # --------------------------------------------------------------------------------------------------------------------- -def datum(h, start): - if not isinstance(h, int): - h = int(h) +def datum(leadtime, start, datetime_obj): + if not isinstance(leadtime, int): + leadtime = int(leadtime) if not isinstance(start, int): start = int(start) - today = datetime.date.today() - today = today.timetuple() - x = datetime.datetime(today.tm_year, today.tm_mon, today.tm_mday, start) - # string = "Forecasttime " + x.strftime("%d.%m. %H:%M") + " UTC " - string1 = "Run: " + x.strftime("%d.%m. %H UTC") - x = datetime.datetime(today.tm_year, today.tm_mon, today.tm_mday, start) + datetime.timedelta(hours=h) - string2 = x.strftime("%A, %d.%m. %H UTC") + today = datetime_obj.timetuple() + x = datetime(today.tm_year, today.tm_mon, today.tm_mday, start) + modelrun_string = "Run: " + x.strftime("%d.%m. %H UTC") + x = datetime(today.tm_year, today.tm_mon, today.tm_mday, start) + timedelta(hours=leadtime) + valid_string = x.strftime("%A, %d.%m. %H UTC") - return string1, string2 + return modelrun_string, valid_string # --------------------------------------------------------------------------------------------------------------------- @@ -58,7 +56,7 @@ def download_nwp(fieldname, datum="20240227", run="00", fp=0, store_path="./"): # --------------------------------------------------------------------------------------------------------------------- -def open_gribfile_single(fieldname, datetime_obj, run, fp, path="./iconnest/"): +def open_gribfile_single(fieldname, datetime_obj, run, fp, path="./modeldata/"): date_string = datetime_obj.strftime("%Y%m%d") nwp_singlelevel = "icon-eu_europe_regular-lat-lon_single-level" @@ -80,9 +78,9 @@ def open_gribfile_single(fieldname, datetime_obj, run, fp, path="./iconnest/"): return data, lats, lons -def open_gribfile_multi(fieldname, lvl, datetime_obj, run, fp, path="./iconnest/"): +def open_icon_gribfile_preslvl(fieldname, lvl, datetime_obj, run, fp, path="./modeldata/"): date_string = datetime_obj.strftime("%Y%m%d") - nwp_modellevel = "icon-eu_europe_regular-lat-lon_model-level" + nwp_modellevel = "icon-eu_europe_regular-lat-lon_pressure-level" # Open the GRIB file grbs = pygrib.open(f"{path}{nwp_modellevel}_{date_string}{run:02d}_{fp:03d}_{lvl}_{fieldname.upper()}.grib2") @@ -100,6 +98,48 @@ def open_gribfile_multi(fieldname, lvl, datetime_obj, run, fp, path="./iconnest/ print("Maximum:", np.nanmax(data)) return data + +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()) + u_fld = np.full(shape, np.nan) + v_fld = np.full(shape, np.nan) + + shape = (model_obj.getnlat(), model_obj.getnlon()) + cape_fld = np.full(shape, np.nan) + 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) + + for par in model_obj.getParamter(): + try: + grb_message = gribidx.select(shortName=par[0], typeOfLevel=par[1], level=par[2])[0] + if par[0] == "cape": + cape_fld = grb_message.values + lats, lons = grb_message.latlons() + else: + idx = pres_levels.index(par[2]) + if par[0] == 'u': + u_fld[idx, :, :] = grb_message.values + elif par[0] == 'v': + v_fld[idx, :, :] = grb_message.values + else: + raise ValueError(f"Unknown Parameter: {par[0]}") + except ValueError as e: + print(f"Error: {e}\t shortName: {par[0]} typeOfLevel: {par[1]} level: {par[2]}") + pass + + return cape_fld, u_fld, v_fld, lats, lons + # --------------------------------------------------------------------------------------------------------------------- diff --git a/test/.coverage b/test/.coverage new file mode 100644 index 0000000..9939b10 Binary files /dev/null and b/test/.coverage differ diff --git a/test/nose2.cfg b/test/nose2.cfg new file mode 100644 index 0000000..84743c4 --- /dev/null +++ b/test/nose2.cfg @@ -0,0 +1,2 @@ +[coverage] +always-on = True diff --git a/test/test_modelinfo.py b/test/test_modelinfo.py new file mode 100755 index 0000000..4cef603 --- /dev/null +++ b/test/test_modelinfo.py @@ -0,0 +1,57 @@ +#!/usr/bin/python3 +import sys +import unittest +from datetime import datetime +# project moduls +sys.path.append('../src/') +from modelinfolib import MODELIFNO + + +class TestMODELIFNO(unittest.TestCase): + def setUp(self): + self.icon_nest = MODELIFNO("ICON Nest", 1377, 657, 0.0625, "pres") + + def test_setrun(self): + run = "15" + self.icon_nest.setrun(run) + self.assertEqual(self.icon_nest.getrun(), int(run)) + + def test_setrundate_with_date_object(self): + test_date = "2024-04-21" + self.icon_nest.setrundate(test_date) + self.assertEqual(self.icon_nest.getrundate(), datetime(2024, 4, 21)) + self.assertEqual(self.icon_nest.getrundate_as_str(fmt="%Y-%m-%d"), "2024-04-21") + + def test_getname(self): + self.assertEqual(self.icon_nest.getname(), "ICON Nest") + + def test_getParamter(self): + self.assertIsNone(self.icon_nest.getParamter()) + + def test_getpoints(self): + self.assertEqual(self.icon_nest.getpoints(), 904689) + + def test_getnlon(self): + self.assertEqual(self.icon_nest.getnlon(), 1377) + + def test_getnlat(self): + self.assertEqual(self.icon_nest.getnlat(), 657) + + def test_getlevels(self): + self.assertEqual(self.icon_nest.getlevels(), [1000, 950, 925, 900, 875, 850, 825, 800, 775, 700, 600, 500, 400, 300]) + + def test_getnlev(self): + self.assertEqual(self.icon_nest.getnlev(), 14) + + def test_getlevtyp(self): + self.assertEqual(self.icon_nest.getlevtyp(), 0.0625) + + def test_getd_grad(self): + self.assertEqual(self.icon_nest.getd_grad(), 0.0625) + + def test_create_plottitle(self): + self.assertEqual(self.icon_nest.create_plottitle(), "Hodographmap of ICON Nest") + + +if __name__ == '__main__': + unittest.main()