Merge pull request #127 from rom-py/125-schism-winds-sflux-generation…

…-produces-misinterpreted-netcdf

125 schism winds sflux generation produces misinterpreted netcdf

notebooks/schism/demo.yaml

@@ -3,31 +3,37 @@ period:
   start: 20230101T00
   end: 20230101T12
   interval: 3600
-run_id: test_schism
+run_id: test_schismcsiro
 config:
   model_type: schismcsiro
   mesbf: 1
   fricc: 0.067
+  param_iof_hydro1: 1 # elevation
+  param_iof_hydro2: 1 # mslp
+  param_iof_hydro14: 1 # wind speed
+  param_iof_hydro16: 1 # surface velocities
+  wwm18: 1 # peak wave direction
+  wwm1: 1 # significant wave height
+  wwm9: 1 # peak period
   grid:
     grid_type: schism
     hgrid:
       id: hgrid
       model_type: data_blob
-      source: ../../tests/schism/test_data/hgrid.gr3
-        # source: ../../tests/schism/test_data/hgrid_20kmto60km_rompyschism_testing.gr3
+       #source: ../../tests/schism/test_data/hgrid.gr3
+      source: ../../tests/schism/test_data/hgrid_20kmto60km_rompyschism_testing.gr3
     manning: 1
   data:
     data_type: schism
     atmos:
       air_1:
         data_type: sflux_air
         source: 
-          model_type: source_file
-          uri: "../../tests/schism/test_data/atmos.nc"
-          model_type: open_dataset
+          model_type: file
+          uri: "../../tests/schism/test_data/era5.nc"
         uwind_name: u10
         vwind_name: v10
-        prmsl_name: mslp
+        prmsl_name: msl
         filter:
           sort: {coords: [latitude]}
         buffer: 5
@@ -41,7 +47,7 @@ config:
           z: depth
         source: 
           uri:  ../../tests/schism/test_data/hycom.nc
-          model_type: open_dataset
+          model_type: file
         variable: surf_el
     tides:
       constituents:

pyproject.toml

@@ -122,3 +122,6 @@ log_cli = true
 log_cli_level = "INFO"
 log_cli_format = "%(asctime)s [%(levelname)8s] %(message)s (%(filename)s:%(lineno)s)"
 log_cli_date_format = "%Y-%m-%d %H:%M:%S"
+
+[tool.black]
+line-length = 88

rompy/schism/data.py

@@ -118,6 +118,18 @@ def ds(self):
         }
         ds.time.encoding["units"] = unit
         ds.time.encoding["calendar"] = "proleptic_gregorian"
+        # open bad dataset
+
+        # SCHISM doesn't like scale_factor and add_offset attributes and requires Float64 values
+        for var in ds.data_vars:
+            # If the variable has scale_factor or add_offset attributes, remove them
+            if "scale_factor" in ds[var].encoding:
+                del ds[var].encoding["scale_factor"]
+            if "add_offset" in ds[var].encoding:
+                del ds[var].encoding["add_offset"]
+            # set the data variable encoding to Float64
+            ds[var].encoding["dtype"] = np.dtypes.Float64DType()
+
         return ds
 
 
@@ -288,6 +300,15 @@ def check_weights(v):
                     f"Relative weights for {variable} do not add to 1.0: {weight}"
                 )
             return v
+        # SCHISM doesn't like scale_factor and add_offset attributes and requires Float64 values
+        for var in ds.data_vars:
+            # If the variable has scale_factor or add_offset attributes, remove them
+            if "scale_factor" in ds[var].encoding:
+                del ds[var].encoding["scale_factor"]
+            if "add_offset" in ds[var].encoding:
+                del ds[var].encoding["add_offset"]
+            # set the data variable encoding to Float64
+            ds[var].encoding["dtype"] = np.dtypes.Float64DType()
 
 
 class SCHISMDataWave(BoundaryWaveStation):
@@ -341,6 +362,20 @@ def get(
         ds.spec.to_ww3(outfile)
         return outfile
 
+    @property
+    def ds(self):
+        """Return the filtered xarray dataset instance."""
+        ds = super().ds
+        for var in ds.data_vars:
+            # If the variable has scale_factor or add_offset attributes, remove them
+            if "scale_factor" in ds[var].encoding:
+                del ds[var].encoding["scale_factor"]
+            if "add_offset" in ds[var].encoding:
+                del ds[var].encoding["add_offset"]
+            # set the data variable encoding to Float64
+            ds[var].encoding["dtype"] = np.dtypes.Float64DType()
+        return ds
+
     def __str__(self):
         return f"SCHISMDataWave"
 
@@ -479,6 +514,15 @@ def boundary_ds(self, grid: SCHISMGrid, time: Optional[TimeRange]) -> xr.Dataset
         if schism_ds.time_series.isnull().any():
             msg = "Some values are null. This will cause SCHISM to crash. Please check your data."
             logger.warning(msg)
+
+        # If the variable has scale_factor or add_offset attributes, remove them
+        # and set the data variable encoding to Float64
+        for var in schism_ds.data_vars:
+            if "scale_factor" in schism_ds[var].encoding:
+                del schism_ds[var].encoding["scale_factor"]
+            if "add_offset" in schism_ds[var].encoding:
+                del schism_ds[var].encoding["add_offset"]
+            schism_ds[var].encoding["dtype"] = np.dtypes.Float64DType()
         return schism_ds
 
 

rompy/schism/utils.py

@@ -1,5 +1,6 @@
 import logging
 import re
+import sys
 from pathlib import Path
 
 import cartopy.crs as ccrs
@@ -81,12 +82,16 @@ def schism_plot(
         contours = np.asarray(contours)
     if varname == "depth" or varname == "z":
         var = z
+    if varname == "wind_speed":
+        var = np.sqrt(schout.wind_speed[:, 0] ** 2 + schout.wind_speed[:, 1] ** 2)
+    if varname == "dahv":
+        var = np.sqrt(schout.dahv[:, 0] ** 2 + schout.dahv[:, 1] ** 2)
     else:
         var = schout[varname]
 
     if len(varscale) == 0:
-        vmin = var.min()
-        vmax = var.max()
+        vmin = var.values.min()
+        vmax = var.values.max()
     else:
         vmin, vmax = varscale
     if project:
@@ -135,6 +140,8 @@ def schism_plot(
     if vectors:
         if re.search("WWM", varname):
             vtype = "waves"
+        if re.search("wind", varname):
+            vtype = "wind"
         else:
             vtype = "currents"
         LonI, LatI, UI, VI = schism_calculate_vectors(ax, schout, vtype=vtype)
@@ -210,11 +217,13 @@ def schism_calculate_vectors(ax, schout, vtype="waves", dX="auto", mask=True):
     elif vtype == "elev" or re.search("curr", vtype):
         idx = np.sqrt(schout.dahv[:, 0] ** 2 + schout.dahv[:, 1] ** 2) > 0.01
         u = schout.dahv[idx, 0]
-        v = schout.dahv[
-            idx, 1
-        ]  # dahv has u and v components, so use index of 1 for v and index of 0 for u
+        v = schout.dahv[idx, 1]
+    elif vtype == "wind":
+        idx = np.ones_like(schout.wind_speed[:, 0], dtype=bool)
+        u = schout.wind_speed[idx, 0]
+        v = schout.wind_speed[idx, 1]
     else:
-        print("*** Warning input vecter data not understood")
+        raise ValueError("*** Warning input vecter data not understood")
     x, y = pUTM55(
         schout.SCHISM_hgrid_node_x.values[idx], schout.SCHISM_hgrid_node_y.values[idx]
     )
@@ -246,26 +255,29 @@ def schism_calculate_vectors(ax, schout, vtype="waves", dX="auto", mask=True):
 
 if __name__ == "__main__":
     # load schism files
-    schfile = "../../notebooks/schism/schism_procedural/test_schism/outputs/schout_1.nc"
+    # schfile = "../../notebooks/schism/schism_procedural/test_schism/outputs/schout_1.nc"
+    # take schfile off the command line
+    schfile = sys.argv[1]
     schout, meshtri = schism_load(schfile)
     lons = schout.SCHISM_hgrid_node_y.values
     lats = schout.SCHISM_hgrid_node_x.values
     # plot gridded fields - elevation
-    for ix, time in enumerate(schout.time[5:8].values):
-        fig, ax = schism_plot(
-            schout,
-            meshtri,
-            "elev",
-            bbox=[145, -25, 155, -16],
-            project=True,
-            plotmesh=True,
-            time=time,
-            mask=False,
-            vectors=True,
-            varscale=(-9, 9),
-            contours=[0],
-        )
-        ax.tick_params(axis="both", which="major", labelsize=24)
-        plt.savefig(f"schism_plot_{ix}.png", dpi=300)
-        plt.show()
-        # plt.close()
+    # for variable in ["elev", "wind_speed", "WWM_1"]:
+    for variable in ["dahv"]:
+        for ix, time in enumerate(schout.time.values):
+            fig, ax = schism_plot(
+                schout,
+                meshtri,
+                variable,
+                bbox=[145, -25, 155, -16],
+                project=True,
+                plotmesh=True,
+                time=time,
+                mask=False,
+                vectors=True,
+                # varscale=(-9, 9),
+                contours=[0],
+            )
+            plt.savefig(f"schism_plot_{variable}_{ix}.png", dpi=300)
+    # plt.show()
+    # plt.close()

tests/schism/test_data/download_era5_datamesh.py

@@ -0,0 +1,59 @@
+# If you are working in a notebook or ipython environment, you may need to uncomment the line below to install the oceanum python library
+#!pip install oceanum
+
+import xarray as xr
+# Define your datamesh token in evironment variables as DATAMESH_TOKEN or insert into argument below.
+from oceanum.datamesh import Connector
+
+datamesh = Connector()
+
+from oceanum.datamesh import Connector
+
+geom = [
+    140,
+    -26,
+    157,
+    -15,
+]
+
+times = ["2023-01-01T00:00:00.000Z", "2023-01-02T00:00:00.000Z"]
+
+# Put your datamesh token in the Jupyterlab settings, or as argument in the constructor below
+datamesh = Connector()
+
+era5_mslp = datamesh.query(
+    {
+        "id": "00bcb40ecbeea137e7350ae9f9f9a3b3",
+        "label": "era5_mslp_00bcb4",
+        "geofilter": {
+            "geom": geom,
+            "type": "bbox",
+            "interp": "linear",
+        },
+        "variables": ["msl"],
+        "datasource": "era5_mslp",
+        "timefilter": {
+            "times": ["2023-01-01T00:00:00.000Z", "2023-01-04T00:00:00.000Z"]
+        },
+        "description": "ECMWF ERA5 global mean sea leve pressure hindcast",
+    }
+)
+era5_wind10m = datamesh.query(
+    {
+        "id": "bdaab8e6341f95e35764bcd1a9f737b5",
+        "label": "era5_wind10m_bdaab8",
+        "geofilter": {
+            "geom": geom,
+            "type": "bbox",
+            "interp": "linear",
+        },
+        "variables": ["u10", "v10"],
+        "datasource": "era5_wind10m",
+        "timefilter": {
+            "times": ["2023-01-01T00:00:00.000Z", "2023-01-04T00:00:00.000Z"]
+        },
+        "description": "ECMWF ERA5 global 10m wind reanalysis",
+    }
+)
+
+xr.merge([era5_mslp, era5_wind10m]).to_netcdf("era5.nc")