Skip to content

Commit

Permalink
GFS and IFS as model
Browse files Browse the repository at this point in the history
  • Loading branch information
@mammatus95
mammatus95 committed Apr 10, 2024
1 parent 1d7efc0 commit 4677c1f
Show file tree
Hide file tree
Showing 6 changed files with 178 additions and 47 deletions.
5 changes: 3 additions & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,9 @@ bash download_script.sh
conda activate HodographMaps

# Plot Hodograph
python3 main.py Basic

python3 main.py IFS
python3 main.py GFS
python3 main.py ICON
cd images
```

Expand Down
9 changes: 6 additions & 3 deletions src/download_script.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,12 +59,15 @@ do
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
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.goessimpgrb2.0p25.f${T}
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

Expand Down
61 changes: 34 additions & 27 deletions src/main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,38 +18,45 @@
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="./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 == "Basic":
cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, 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
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()
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, 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"])
else:
print("Wrong command line argument")
exit(-1)
if program_mode == "Test":
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 == "Basic":
cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, 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_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, 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"])
else:
print("Wrong command line argument")
exit(-1)

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

Expand Down
91 changes: 84 additions & 7 deletions src/modelinfolib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,15 @@
d_grad = 0.0625
# IFS
points = 405900
nlon = 900
nlat = 451
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
Expand All @@ -34,6 +34,78 @@
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),
('u', 'isobaricInhPa', 250),
('u', 'isobaricInhPa', 200),
('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),
('v', 'isobaricInhPa', 250),
('v', 'isobaricInhPa', 200),
('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),
('u', 'isobaricInhPa', 250),
('u', 'isobaricInhPa', 200),
('v', 'isobaricInhPa', 1000),
('v', 'isobaricInhPa', 925),
('v', 'isobaricInhPa', 850),
('v', 'isobaricInhPa', 700),
('v', 'isobaricInhPa', 600),
('v', 'isobaricInhPa', 500),
('v', 'isobaricInhPa', 400),
('v', 'isobaricInhPa', 300),
('v', 'isobaricInhPa', 250),
('v', 'isobaricInhPa', 200),
('cape', 'entireAtmosphere', 0),
# ('10u', 'heightAboveGround', 10),
# ('10v', 'heightAboveGround', 10),
# ('100u', 'heightAboveGround', 100),
# ('100v', 'heightAboveGround', 100),
]

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


class MODELIFNO:
Expand All @@ -45,16 +117,18 @@ def __init__(self, modelname, nlon, nlat, d_grad, levtyp):
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.parlist = None
elif modelname == "IFS":
self.levels = [1000, 925, 850, 700, 600, 500, 400, 300, 250, 200]
self.parlist = par_list_ifs
elif modelname == "GFS":
self.levels = [1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 500, 400, 300, 250, 200]
self.parlist = par_list_gfs
else:
self.levels = [1000, 850, 700, 600, 500, 400, 300, 250, 200]
self.d_grad = d_grad
self.levtyp = levtyp


def __str__(self):
return (
f"Model Information: {self.modelname}\nPoints: {self.points}\n"
Expand All @@ -66,6 +140,9 @@ def __str__(self):
def getname(self):
return self.modelname

def getParamter(self):
return self.parlist

def getpoints(self):
return self.points

Expand Down Expand Up @@ -96,5 +173,5 @@ def getd_grad(self):

# Example usage:
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")
ifs = MODELIFNO("IFS", 1440, 721, 0.25, "pres")
gfs = MODELIFNO("GFS", 1440, 721, 0.25, "pres")
14 changes: 8 additions & 6 deletions src/plotlib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -235,9 +235,11 @@ def basic_plot(cape_fld, u, v, pres_levels, lats, lons, hour, start, titel='CAPE
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):
if np.mean(cape_fld[i-1:i+1, j-1:j+1]) > threshold:
#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:
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()
Expand All @@ -249,9 +251,9 @@ def basic_plot(cape_fld, u, v, pres_levels, lats, lons, hour, start, titel='CAPE
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('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_ce_{hour}.{imfmt}"
Expand Down
45 changes: 43 additions & 2 deletions src/utilitylib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,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"

Expand All @@ -80,7 +80,7 @@ def open_gribfile_single(fieldname, datetime_obj, run, fp, path="./iconnest/"):
return data, lats, lons


def open_icon_gribfile_preslvl(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_pressure-level"

Expand All @@ -100,6 +100,47 @@ def open_icon_gribfile_preslvl(fieldname, lvl, datetime_obj, run, fp, path="./ic
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")

# 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()

# 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] # takes the matching grib message
print(grb_message)
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

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


Expand Down

0 comments on commit 4677c1f

Please sign in to comment.