Skip to content

Commit

Permalink
Lint (#1)
Browse files Browse the repository at this point in the history
* lint 1

* lint 2

* lint 3

* Create LICENSE

* Add License to README.md
  • Loading branch information
mammatus95 authored Apr 9, 2024
1 parent eb56a2a commit 977ab27
Show file tree
Hide file tree
Showing 10 changed files with 212 additions and 158 deletions.
21 changes: 21 additions & 0 deletions LICENSE
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
MIT License

Copyright (c) 2024 Morten Kretschmer

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
4 changes: 4 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -42,3 +42,7 @@ python3 main.py Sounding
## Datasource
- ICON Nest: https://opendata.dwd.de/weather/nwp/icon-eu/
- IFS: https://www.ecmwf.int/en/forecasts/datasets/open-data

## License

This project is licensed under the terms of the MIT license. See the [LICENSE](LICENSE) file for details.
Binary file modified images/hodographmap_ce_9.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
22 changes: 14 additions & 8 deletions src/creategridinfo.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,23 +3,29 @@

import numpy as np

path_icon='./'
path_icon = './'

if (os.getcwd() != path_icon):
os.chdir(path_icon)


def creategridinfo(d_grad=0.0625, filename="gitter_icon"):
lonmin = -23.5
lonmax = 45.0
latmin = 29.5
latmax = 70.5
xsize=np.size(np.arange(lonmin, lonmax+d_grad, d_grad))
ysize=np.size(np.arange(latmin, latmax+d_grad, d_grad))
test_string='#\n# gridID 1\n#\ngridtype = lonlat\ngridsize = %d\nxname = lon\nxlongname = longitude\nxunits = degrees_east\nyname = lat\nylongname = latitude\nyunits = degrees_north\nxsize = %d\nysize = %d\nxfirst = -23.5\nxinc = %.1f\nyfirst = 29.5\nyinc = %.1f' %(xsize*ysize, xsize, ysize, d_grad, d_grad)
fd = open(filename, "x")
lonmax = 45.0
latmin = 29.5
latmax = 70.5
xsize = np.size(np.arange(lonmin, lonmax+d_grad, d_grad))
ysize = np.size(np.arange(latmin, latmax+d_grad, d_grad))
test_string = (f"#\n# gridID 1\n#\ngridtype = lonlat\ngridsize = {xsize*ysize}\nxname = lon\nxlongname = longitude"
f"\nxunits = degrees_east\nyname = lat\nylongname = latitude\nyunits = degrees_north"
f"\nxsize = {xsize}\nysize = {ysize}\nxfirst = -23.5\nxinc = {d_grad}\n"
f"yfirst = 29.5\nyinc = {d_grad}"
)
fd = open(filename, "x")
fd.write(test_string)
fd.close()


creategridinfo()
creategridinfo(d_grad=0.12, filename="gitter_icon1.2")
creategridinfo(d_grad=0.2, filename="gitter_icon2")
189 changes: 102 additions & 87 deletions src/main.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,73 +12,22 @@
import plotlib
import modelinfolib as model
model = model.icon_nest
# ---------------------------------------------------------------------------------------------------------------------

def main():
config = ut.load_yaml('config.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')

# Add optional argument
parser.add_argument('-f', '--field', type=str,
help='Which background field: CAPE ML, CAPE CON or LPI ')

parser.add_argument('-d', '--date', type=str,
help='Date Format YYYY-MM-DD')

parser.add_argument('-fp', '--fp', type=str,
help='Leadtime or forecast periode')

parser.add_argument('-r', '--run', type=str,
help='Run (0,6,12,18,.. .etc)')

# Parse the command line arguments
args = parser.parse_args()

if args.field == None:
args.field = "CAPE ML"

if (args.field != "CAPE ML") and (args.field != "CAPE CON") and (args.field != "LPI") and (args.field != "WMAXSHEAR"):
print("Unknown field")
exit(-1)

if args.date == None:
rundate = datetime.strptime(config["default_date"], "%Y-%m-%d")
else:
rundate = datetime.strptime(args.date, "%Y-%m-%d")

if args.fp == None:
fp = config["fp"]
else:
fp = args.fp

if args.run == None:
run = config["run"]
else:
run = args.run

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

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

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

#ut.download_nwp(fieldname, datum="20240227", run="00", fp=0, store_path="./")
def run(program_mode, fieldname, rundate, model_run, fp):
config = ut.load_yaml('config.yml')


if args.mode == "Test":
cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, run, fp, path="./iconnest/")
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, run, titel='CAPE')
plotlib.test_plot(cape_fld, lats, lons, fp, model_run, titel='CAPE')

elif args.mode == "Sounding":
cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, run, fp, path="./iconnest/")
elif program_mode == "Sounding":
cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./iconnest/")

steps=config["steps"]
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)
Expand All @@ -95,20 +44,20 @@ def main():

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, run, fp, path="./iconnest/")
q_fld[lvl_idx,:,:] = ut.open_gribfile_multi("QV", level, rundate, run, fp, path="./iconnest/")
p_fld[lvl_idx,:,:] = ut.open_gribfile_multi("P", level, rundate, run, fp, path="./iconnest/")
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, run, titel='CAPE')
elif args.mode == "Basic":
cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, run, fp, path="./iconnest/")
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"]
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)
Expand All @@ -123,20 +72,22 @@ def main():

lvl_idx = 0
for level in range(config["levels"][0], config["levels"][1], steps):
u_fld[lvl_idx,:,:] = ut.open_gribfile_multi("U", level, rundate, run, fp, path="./iconnest/")
v_fld[lvl_idx,:,:] = ut.open_gribfile_multi("V", level, rundate, run, fp, path="./iconnest/")
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

print(np.nanmean(u_fld, axis=(1,2)))
plotlib.basic_plot (cape_fld, u_fld, v_fld, lats, lons, fp, run, titel='CAPE', threshold=config["threshold"])
plotlib.basic_plot_custarea (cape_fld, u_fld, v_fld, lats, lons, fp, run, titel='CAPE', threshold=config["threshold"])
elif args.mode == "Nixon":
cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, run, fp, path="./iconnest/")

steps=config["steps"]

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)
Expand All @@ -155,26 +106,90 @@ def main():

lvl_idx = 0
for level in range(config["levels"][0], config["levels"][1], steps):
u_fld[lvl_idx,:,:] = ut.open_gribfile_multi("U", level, rundate, run, fp, path="./iconnest/")
v_fld[lvl_idx,:,:] = ut.open_gribfile_multi("V", level, rundate, run, fp, path="./iconnest/")
h_fld[lvl_idx,:,:] = ut.open_gribfile_multi("H", level, rundate, run, fp, path="./iconnest/")
p_fld[lvl_idx,:,:] = ut.open_gribfile_multi("P", level, rundate, run, fp, path="./iconnest/")
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)))
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,:,:])
du = np.subtract(u_fld[30, :, :], u_fld[0, :, :])
dv = np.subtract(v_fld[30, :, :], v_fld[0, :, :])
dls_fld = np.sqrt(np.add(np.square(du), np.square(dv)))
plotlib.nixon_proj(cape_fld, dls_fld, u_fld, v_fld, p_fld, h_fld, lats, lons, fp, run, imfmt="png")
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)

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


def main():
config = ut.load_yaml('config.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')

# Add optional argument
parser.add_argument('-f', '--field', type=str,
help='Which background field: CAPE ML, CAPE CON or LPI ')

parser.add_argument('-d', '--date', type=str,
help='Date Format YYYY-MM-DD')

parser.add_argument('-fp', '--fp', type=str,
help='Leadtime or forecast periode')

parser.add_argument('-r', '--run', type=str,
help='Run (0,6,12,18,.. .etc)')

# Parse the command line arguments
args = parser.parse_args()

if args.field is None:
args.field = "CAPE ML"

if (args.field != "CAPE ML") and (args.field != "CAPE CON") and (args.field != "LPI") and (args.field != "WMAXSHEAR"):
print("Unknown field")
exit(-1)

if args.date is None:
rundate = datetime.strptime(config["default_date"], "%Y-%m-%d")
else:
rundate = datetime.strptime(args.date, "%Y-%m-%d")

if args.fp is None:
fp = config["fp"]
else:
fp = args.fp

if args.run is None:
model_run = config["run"]
else:
model_run = args.run

if args.mode != "Test" or args.mode != "Sounding" or args.mode != "Basic" or args.mode != "Nixon":
print("Unknown Mode. Exit program.")
exit(0)
else:
program_mode = args.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)

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


if __name__ == "__main__":
main()

28 changes: 17 additions & 11 deletions src/meteolib.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#!/usr/bin/python3
################################
#
# Meteorology library
# Meteorology library
#
################################

Expand All @@ -15,11 +15,12 @@
GRAV = 9.80665 # Gravity
TOL = 1e-10 # Floating Point Tolerance

def uv2spddir(u,v):

direction = np.rad2deg(np.arctan2(-u,v))
def uv2spddir(u, v):

if type(direction) == np.ndarray:
direction = np.rad2deg(np.arctan2(-u, v))

if isinstance(direction, np.ndarray):
direction[np.where(direction < 0)] += 360
else:
if direction < 0:
Expand All @@ -39,6 +40,7 @@ def q_to_mixrat(q):
"""
return q/(1.0-q)


def temp_at_mixrat(w, p):
"""
Returns the temperature in K of air at the given mixing ratio in g/kg and pressure in hPa
Expand All @@ -55,11 +57,15 @@ def temp_at_mixrat(w, p):
"""

# Constants Used
c1 = 0.0498646455 ; c2 = 2.4082965 ; c3 = 7.07475
c4 = 38.9114 ; c5 = 0.0915 ; c6 = 1.2035
c1 = 0.0498646455
c2 = 2.4082965
c3 = 7.07475
c4 = 38.9114
c5 = 0.0915
c6 = 1.2035

x = np.log10(w * p / (622. + w))
x = (np.power(10.,((c1 * x) + c2)) - c3 + (c4 * np.power((np.power(10,(c5 * x)) - c6),2)))
x = (np.power(10., ((c1 * x) + c2)) - c3 + (c4 * np.power((np.power(10, (c5 * x)) - c6), 2)))
return x


Expand Down Expand Up @@ -92,20 +98,20 @@ def non_parcel_bunkers_motion_experimental(u, v, ps, i_500m, i_5km, i_6km):
"""
d=7.5
d = 7.5

# sfc-500m Mean Wind
mnu500m, mnv500m = mean_wind(u[:i_500m], v[:i_500m], ps[:i_500m])

## 5.5km-6.0km Mean Wind
# 5.5km-6.0km Mean Wind
mnu5500m_6000m, mnv5500m_6000m = mean_wind(u[i_5km:i_6km], v[i_5km:i_6km], ps[i_5km:i_6km])

# shear vector of the two mean winds
shru = mnu5500m_6000m - mnu500m
shrv = mnv5500m_6000m - mnv500m

# SFC-6km Mean Wind
mnu6, mnv6 = mean_wind(u[:i_6km], v[:i_6km], ps[:i_6km])
mnu6, mnv6 = mean_wind(u[:i_6km], v[:i_6km], ps[:i_6km])

# Bunkers Right Motion
tmp = d / np.sqrt(shru*shru + shrv*shrv)
Expand All @@ -114,4 +120,4 @@ def non_parcel_bunkers_motion_experimental(u, v, ps, i_500m, i_5km, i_6km):
lstu = mnu6 - (tmp * shrv)
lstv = mnv6 + (tmp * shru)

return rstu, rstv, lstu, lstv, mnu6, mnv6
return rstu, rstv, lstu, lstv, mnu6, mnv6
Loading

0 comments on commit 977ab27

Please sign in to comment.