Code improvements

.gitignore

@@ -1,7 +1,9 @@
 # ignore ALL .log files
 
-.ipynb_checkpoints/
-*/.ipynb_checkpoints/
+**__pycache__
+
+**.ipynb_checkpoints/
+
 
 */*.png
 */*/*.png

README.md

@@ -16,9 +16,9 @@ You need to addg Conda initialization to your etc/profile, as well.
 
 ![Example](images/example.png)
 
-![Map of Hodographs](images/hodographmap_area.png)
+![Map of Hodographs IFS](images/hodographmap_IFS.png)
 
-![Map of Soundings](images/soundingmap_example.png)
+![Map of Hodographs GFS](images/hodographmap_GFS.png)
 
 
 ## How to run
@@ -35,16 +35,16 @@ python3 main.py ICON
 cd images
 ```
 
+## Datasource
+- [ICON Nest](https://opendata.dwd.de/weather/nwp/icon-eu/)
+- [IFS](https://www.ecmwf.int/en/forecasts/datasets/open-data)
+- [GFS (NOAA)](https://www.nco.ncep.noaa.gov/pmb/products/gfs/)
+
 ## Cartopy?
 
 - https://github.com/mammatus95/cartopy
 - https://scitools.org.uk/cartopy/docs/latest/#
 
-## Datasource
-- ICON Nest: https://opendata.dwd.de/weather/nwp/icon-eu/
-- IFS: https://www.ecmwf.int/en/forecasts/datasets/open-data
-- GFS: https://www.nco.ncep.noaa.gov/pmb/products/gfs/
-
 ## License
 
 This project is licensed under the terms of the MIT license. See the [LICENSE](LICENSE) file for details.

src/main.py

@@ -15,23 +15,29 @@
 # ---------------------------------------------------------------------------------------------------------------------
 
 
-def run(model_obj, program_mode, fieldname, rundate, model_run, fp):
+def run(model_obj, fp):
     config = ut.load_yaml('config.yml')
-    print(model_obj)
+
     if model_obj.getname() == "IFS" or model_obj.getname() == "GFS":
-        cape_fld, u_fld, v_fld, lats, lons = ut.open_gribfile_preslvl(model_obj, rundate, model_run, fp, path="./modeldata/")
-        pres_levels = model_obj.getlevels()
+        cape_fld, u_fld, v_fld, lats, lons = ut.open_gribfile_preslvl(model_obj, fp, path="./modeldata/")
+
         print(np.nanmean(u_fld, axis=(1, 2)))
-        plotlib.basic_plot(cape_fld, u_fld, v_fld, pres_levels, lats, lons, fp, model_run,
-                           titel='CAPE', threshold=config["threshold"])
+        plotlib.basic_plot(model_obj, cape_fld, u_fld, v_fld, lats, lons, fp,
+                           threshold=config["threshold"])
     else:
+        # program_mode
+        program_mode = config["program_mode"]
+
+        # replace space with underscores
+        fieldname = "CAPE_ML"  #args.field.replace(" ", "_")
+        rundate = model_obj.getrundate()
         if program_mode == "Test":
-            cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./modeldata/")
+            cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_obj.getrun(), fp, path="./modeldata/")
             assert cape_fld.shape == (model_obj.getnlat(), model_obj.getnlon()), "Shape inconsistency"
-            plotlib.test_plot(cape_fld, lats, lons, fp, model_run, titel='CAPE')
+            plotlib.test_plot(cape_fld, lats, lons, fp, model_obj.getrun(), titel='CAPE')
 
         elif program_mode == "Basic":
-            cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_run, fp, path="./modeldata/")
+            cape_fld, lats, lons = ut.open_gribfile_single(fieldname, rundate, model_obj.getrun(), fp, path="./modeldata/")
 
             nlvl = int(model_obj.getnlev())
             u_fld = np.empty(nlvl*model_obj.getpoints()).reshape((nlvl, model_obj.getnlat(), model_obj.getnlon()))
@@ -42,18 +48,18 @@ def run(model_obj, program_mode, fieldname, rundate, model_run, fp):
             lvl_idx = 0
             pres_levels = model_obj.getlevels()
             for level in pres_levels:
-                u_fld[lvl_idx, :, :] = ut.open_icon_gribfile_preslvl("U", level, rundate, model_run, fp, path="./modeldata/")
-                v_fld[lvl_idx, :, :] = ut.open_icon_gribfile_preslvl("V", level, rundate, model_run, fp, path="./modeldata/")
+                u_fld[lvl_idx, :, :] = ut.open_icon_gribfile_preslvl("U", level, rundate, model_obj.getrun(), fp, path="./modeldata/")
+                v_fld[lvl_idx, :, :] = ut.open_icon_gribfile_preslvl("V", level, rundate, model_obj.getrun(), fp, path="./modeldata/")
 
                 lvl_idx += 1
                 if lvl_idx >= nlvl:
                     break
 
             print(np.nanmean(u_fld, axis=(1, 2)))
-            plotlib.basic_plot(cape_fld, u_fld, v_fld, pres_levels, lats, lons, fp, model_run,
-                               titel='CAPE', threshold=config["threshold"])
-            plotlib.basic_plot_custarea(cape_fld, u_fld, v_fld, pres_levels, lats, lons, fp, model_run,
-                                        titel='CAPE', threshold=config["threshold"])
+            plotlib.basic_plot(model_obj, cape_fld, u_fld, v_fld, lats, lons, fp,
+                               threshold=config["threshold"])
+            plotlib.basic_plot_custarea(model_obj, cape_fld, u_fld, v_fld, lats, lons, fp,
+                                        threshold=config["threshold"])
         else:
             print("Wrong command line argument")
             exit(-1)
@@ -62,7 +68,6 @@ def run(model_obj, program_mode, fieldname, rundate, model_run, fp):
 
 
 def main():
-    config = ut.load_yaml('config.yml')
     run_config = ut.load_yaml('run.yml')
 
     # command line arguments
@@ -94,21 +99,6 @@ def main():
         print("Unknown field")
         exit(-1)
 
-    if args.date is None:
-        rundate = datetime.strptime(run_config["default_date"], "%Y-%m-%d")
-    else:
-        rundate = datetime.strptime(args.date, "%Y-%m-%d")
-
-    if args.fp is None:
-        fp = run_config["fp"]
-    else:
-        fp = args.fp
-
-    if args.run is None:
-        model_run = run_config["run"]
-    else:
-        model_run = args.run
-
     if args.Model is None:
         model_obj = model.MODELIFNO("ICON EU", 1377, 657, 0.0625, "pres")
     elif "ICON" in args.Model:
@@ -121,15 +111,24 @@ def main():
         print("Unkown model! Exit.")
         exit(0)
 
-    # program_mode
-    program_mode = config["program_mode"]
+    if args.date is None:
+        model_obj.setrundate(datetime.strptime(run_config["default_date"], "%Y-%m-%d"))
+    else:
+        model_obj.setrundate(datetime.strptime(args.date, "%Y-%m-%d"))
 
-    # replace space with underscores
-    fieldname = args.field.replace(" ", "_")
+    if args.fp is None:
+        fp = run_config["fp"]
+    else:
+        fp = args.fp
+
+    if args.run is None:
+        model_obj.setrun(run_config["run"])
+    else:
+        model_obj.setrun(args.run)
 
-    print(f"\nDate: {rundate}\n Arguments: {args} \nConfig-File: {config}\n\n")
+    print(f"\nArguments: {args}\n{model_obj}")
 
-    run(model_obj, program_mode, fieldname, rundate, model_run, fp)
+    run(model_obj, fp)
 
 # ---------------------------------------------------------------------------------------------------------------------
 

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

src/modelinfolib.py

@@ -1,4 +1,7 @@
 #!/usr/bin/python3
+from datetime import datetime, date
+from utilitylib import load_yaml
+
 """
 # ICON Nest
 points = 904689
@@ -51,8 +54,6 @@
                ('u', 'isobaricInhPa',  500),
                ('u', 'isobaricInhPa',  400),
                ('u', 'isobaricInhPa',  300),
-               ('u', 'isobaricInhPa',  250),
-               ('u', 'isobaricInhPa',  200),
                ('v', 'isobaricInhPa', 1000),
                ('v', 'isobaricInhPa',  975),
                ('v', 'isobaricInhPa',  950),
@@ -67,8 +68,6 @@
                ('v', 'isobaricInhPa',  500),
                ('v', 'isobaricInhPa',  400),
                ('v', 'isobaricInhPa',  300),
-               ('v', 'isobaricInhPa',  250),
-               ('v', 'isobaricInhPa',  200),
                ('cape', 'pressureFromGroundLayer', 9000),  # 9000 18000 25500
 # ('cape', 'surface', 0),
 # ('cape', 'pressureFromGroundLayer', 9000),  # 9000 18000 25500
@@ -86,8 +85,6 @@
                ('u', 'isobaricInhPa',  500),
                ('u', 'isobaricInhPa',  400),
                ('u', 'isobaricInhPa',  300),
-               ('u', 'isobaricInhPa',  250),
-               ('u', 'isobaricInhPa',  200),
                ('v', 'isobaricInhPa', 1000),
                ('v', 'isobaricInhPa',  925),
                ('v', 'isobaricInhPa',  850),
@@ -96,8 +93,6 @@
                ('v', 'isobaricInhPa',  500),
                ('v', 'isobaricInhPa',  400),
                ('v', 'isobaricInhPa',  300),
-               ('v', 'isobaricInhPa',  250),
-               ('v', 'isobaricInhPa',  200),
                ('cape', 'entireAtmosphere', 0),
 # ('10u',  'heightAboveGround', 10),
 # ('10v', 'heightAboveGround',  10),
@@ -111,32 +106,63 @@
 class MODELIFNO:
 
     def __init__(self, modelname, nlon, nlat, d_grad, levtyp):
+        config = load_yaml('config.yml')
         self.modelname = modelname
         self.points = nlon*nlat
         self.nlon = nlon
         self.nlat = nlat
         if "ICON" in modelname:
-            self.levels = [1000, 950, 925, 900, 875, 850, 825, 800, 775, 700, 600, 500, 400, 300, 250, 200]
+            self.levels = [1000, 950, 925, 900, 875, 850, 825, 800, 775, 700, 600, 500, 400, 300]
             self.parlist = None
+            self.hodo_interval_lat = range(272, 415, 12)
+            self.hodo_interval_lon = range(420, 670, 15)
         elif modelname == "IFS":
-            self.levels = [1000, 925, 850, 700, 600, 500, 400, 300, 250, 200]
+            self.levels = [1000, 925, 850, 700, 600, 500, 400, 300]
             self.parlist = par_list_ifs
+            self.hodo_interval_lat = range(143, 174, 3)
+            self.hodo_interval_lon = range(731, 794, 3)
         elif modelname == "GFS":
-            self.levels = [1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 500, 400, 300, 250, 200]
+            self.levels = [1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 500, 400, 300]
             self.parlist = par_list_gfs
+            self.hodo_interval_lat = range(143, 174, 3)
+            self.hodo_interval_lon = range(11, 74, 3)
         else:
-            self.levels = [1000, 850, 700, 600, 500, 400, 300, 250, 200]
+            self.levels = [1000, 850, 700, 600, 500, 400, 300]
         self.d_grad = d_grad
         self.levtyp = levtyp
+        self.run = -99
+        self.rundate = date.today()
 
     def __str__(self):
         return (
-                f"Model Information: {self.modelname}\nPoints: {self.points}\n"
-                f"Number of Longitudes: {self.nlon}\nNumber of Latitudes: {self.nlat}\n"
+                f"Model Information: {self.modelname}  Run: {self.run} Date: {self.rundate}\n"
+                f"Points: {self.points}\n"
+                f"Number of Longitudes: {self.nlon}\tNumber of Latitudes: {self.nlat}\n"
                 f"Horizontal resolution: {self.d_grad}\n"
                 f"Number of Levels: {len(self.levels)}\tLeveltyps: {self.levtyp}\n"
                )
 
+    def setrun(self, run):
+        if isinstance(run, int):
+            self.run = run
+        else:
+            self.run = int(run)
+
+    def setrundate(self, rundate):
+        if isinstance(rundate, datetime):
+            self.rundate = rundate
+        else:
+            self.rundate = datetime.strptime(rundate, "%Y-%m-%d")
+
+    def getrun(self):
+        return self.run
+
+    def getrundate(self):
+        return self.rundate
+
+    def getrundate_as_str(self, fmt="%Y%m%d"):
+        return self.rundate.strftime(fmt)
+
     def getname(self):
         return self.modelname
 
@@ -170,8 +196,10 @@ def getlevtyp(self):
     def getd_grad(self):
         return self.d_grad
 
+    def create_plottitle(self):
+        return f"Hodographmap of {self.modelname}"
 
 # Example usage:
-icon_nest = MODELIFNO("ICON EU", 1377, 657, 0.0625, "pres")
+icon_nest = MODELIFNO("ICON Nest", 1377, 657, 0.0625, "pres")
 ifs = MODELIFNO("IFS", 1440, 721, 0.25, "pres")
-gfs = MODELIFNO("GFS", 1440, 721, 0.25, "pres")
+gfs = MODELIFNO("GFS", 1440, 721, 0.25, "pres")