Skip to content

Commit

Permalink
Pressure Levels for ICON
Browse files Browse the repository at this point in the history
  • Loading branch information
@mammatus95
mammatus95 committed Apr 9, 2024
1 parent e9dc14b commit 1d7efc0
Show file tree
Hide file tree
Showing 9 changed files with 109 additions and 247 deletions.
3 changes: 1 addition & 2 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,5 +16,4 @@
*/*.pyc
*/*/*.pyc

src/invariant
src/iconnest
src/modeldata
3 changes: 0 additions & 3 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,6 @@ conda activate HodographMaps
# Plot Hodograph
python3 main.py Basic

# Plot Soundings
python3 main.py Sounding

cd images
```

Expand Down
Binary file removed images/soundingmap_example.png
View file
Binary file not shown.
5 changes: 1 addition & 4 deletions src/config.yml
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,6 @@
levels: [74,20]
steps: 1 # level steps
program_mode : Basic # Test Basic

hodomap : True
soundmap : True

threshold: 0

Expand All @@ -16,4 +14,3 @@ customize: # east de
lat1: 50.1
lat2: 54.8

# 10.77, 18.92, 49.80, 55.06
44 changes: 30 additions & 14 deletions src/download_script.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand All @@ -13,21 +13,27 @@ mkdir -p ./images
# select run
R=0
# 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")

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
do
Expand All @@ -37,21 +43,31 @@ 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/0p4-beta/oper/${D}$(printf "%02d" "$R")0000-${X}h-oper-fc.grib2
wget ${ifs_file} -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

# gfs
gfs_file=${gfs_model_pfad}/gfs.${D}/$(printf "%02d" "$R")/atmos/gfs.t$(printf "%02d" "$R")z.goessimpgrb2.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
Expand All @@ -65,4 +81,4 @@ do
done

# remove nwp files
rm -rf ${store_path}
#rm -rf ${store_path}
125 changes: 30 additions & 95 deletions src/main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,116 +11,42 @@
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, program_mode, fieldname, rundate, model_run, fp):
config = ut.load_yaml('config.yml')

print(model_obj)
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"
cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, 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')

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/")
cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./modeldata/")

steps = config["steps"]
nlvl = int(model.getnlev()/steps)
u_fld = np.empty(nlvl*model.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon()))
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.getpoints()).reshape((nlvl, model.getnlat(), model.getnlon()))
v_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.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/")
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/")

lvl_idx += 1
if lvl_idx >= nlvl:
break

print(np.nanmean(u_fld, axis=(1, 2)))
plotlib.basic_plot(cape_fld, u_fld, v_fld, lats, lons, fp, model_run,
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, lats, lons, fp, model_run,
plotlib.basic_plot_custarea(cape_fld, u_fld, v_fld, pres_levels, 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")
else:
print("Wrong command line argument")
exit(-1)
Expand All @@ -135,8 +61,8 @@ def main():
# 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,
Expand Down Expand Up @@ -176,18 +102,27 @@ def main():
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)
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:
program_mode = args.mode
print("Unkown model! Exit.")
exit(0)

# program_mode
program_mode = config["program_mode"]

# replace space with underscores
fieldname = args.field.replace(" ", "_")

print(f"\nDate: {rundate}\n Arguments: {args} \nConfig-File: {config}\n\n")

run(program_mode, fieldname, rundate, model_run, fp)
run(model_obj, program_mode, fieldname, rundate, model_run, fp)

# ---------------------------------------------------------------------------------------------------------------------

Expand Down
49 changes: 31 additions & 18 deletions src/modelinfolib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
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
Expand All @@ -15,7 +15,7 @@
points = 405900
nlon = 900
nlat = 451
nlev =
lev = [1000, 925, 850, 700, 600, 500, 400, 300, 250, 200]
lonmin = 0
lonmax = 359
latmin = -90
Expand All @@ -26,7 +26,7 @@
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
Expand All @@ -38,32 +38,34 @@

class MODELIFNO:

def __init__(self, nlon, nlat, d_grad, nlev, levtyp):
def __init__(self, modelname, nlon, nlat, d_grad, levtyp):
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, 250, 200]
elif modelname == "IFS":
self.levels = [1000, 925, 850, 700, 600, 500, 400, 300, 250, 200]
elif modelname == "GFS":
self.levels = [1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 500, 400, 300, 250, 200]
else:
self.levels = [1000, 850, 700, 600, 500, 400, 300, 250, 200]
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)
"""

def __str__(self):
return (
f"Model Information:\nPoints: {self.points}\n"
f"Model Information: {self.modelname}\nPoints: {self.points}\n"
f"Number of Longitudes: {self.nlon}\nNumber 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 getname(self):
return self.modelname

def getpoints(self):
return self.points

Expand All @@ -73,8 +75,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
Expand All @@ -84,4 +95,6 @@ def getd_grad(self):


# Example usage:
icon_nest = MODELIFNO(1377, 657, 0.0625, 54, "model") # lowest 74 and we download till 20
icon_nest = MODELIFNO("ICON EU", 1377, 657, 0.0625, "pres")
ifs = MODELIFNO("IFS", 451, 900, 0.4, "pres")
gfs = MODELIFNO("GFS", 721, 1440, 0.25, "pres")
Loading

0 comments on commit 1d7efc0

Please sign in to comment.