[run-ex] Added method to identify and select (add standard_name attribute) surface variables (wind) from GRIB datasets. Added example_grib_cf to illustrate use with global NCAR GFS dataset.

Author: knutfrode

examples/example_grib_cf.py

@@ -0,0 +1,36 @@
+#!/usr/bin/env python
+"""
+Making GRIB dataset CF-compatible
+=================================
+"""
+
+from datetime import datetime, timedelta
+from opendrift.models.windblow import WindBlow
+
+#%%
+# Comparing / using two different wind datasets
+#
+# The first dataset (UCAR) is a GRIB dataset served by Thredds as a netCDF dataset, but with GRIB attributes
+#   CF standard_name are added for identified surface variables (wind) with this method:
+#     https://opendrift.github.io/_modules/opendrift/readers.html#add_standard_name_for_surface_grib_variables
+#   This should be used with care, making sure that only surface (10m height) variables are used
+#
+# The second dataset (PACIOOS) is a netCDF4 dataset served by Thredds with CF attributes
+wind_datasets = [
+    'https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p25deg/Best',
+    'https://pae-paha.pacioos.hawaii.edu/thredds/dodsC/ncep_global/NCEP_Global_Atmospheric_Model_best.ncd'
+    ]
+
+simulations = []
+
+time = datetime.now() - timedelta(days=5)
+for wind in wind_datasets:
+    o = WindBlow(loglevel=20)
+    o.add_readers_from_list(wind)
+    o.seed_elements(lon=3, lat=60, time=time, number=10, radius=10000)
+    o.run(duration=timedelta(hours=48))
+    simulations.append(o)
+
+#%%
+# Some differences are expected as datasets may origin from different runs
+simulations[0].plot(compare=simulations[1:], lscale='i', legend=['GRIB converted', 'netCDF CF compatible'])

opendrift/models/windblow.py

@@ -35,7 +35,7 @@ class WindBlow(OpenDriftSimulation):
 
     def __init__(self, *args, **kwargs):
         super(WindBlow, self).__init__(*args, **kwargs)
-        self._set_config_default('drift:max_speed', 12)
+        self._set_config_default('drift:max_speed', 25)
 
     def update(self):
 

opendrift/readers/__init__.py

@@ -25,6 +25,7 @@
 import logging; logger = logging.getLogger(__name__)
 import glob
 import json
+import numpy as np
 import opendrift
 import xarray as xr
 import copernicusmarine
@@ -85,6 +86,40 @@ def open_dataset_opendrift(source, zarr_storage_options=None, open_mfdataset_opt
 
     return ds
 
+def add_standard_name_for_surface_grib_variables(ds):
+    """Try to identify surface variables from GRIB datasets
+       and add standard_name attribute"""
+    mapping = {
+        'wind_speed': {'Grib2_Parameter': np.array([0, 2, 1])},
+        'x_wind': {'Grib2_Parameter': np.array([0, 2, 2])},
+        'y_wind': {'Grib2_Parameter': np.array([0, 2, 3])},
+        }
+
+    for var_name, var in ds.data_vars.items():
+        if any(forbidden.lower() in var_name.lower() for forbidden in ('percentile', 'pctl')):
+            logger.debug(f'Skipping percentile variable {var_name}')
+            continue
+        if 'standard_name' in var.attrs:
+            continue
+        if var.attrs.get('Grib2_Level_Type', None) != 103:
+            continue  # Not surface variable
+        for standard_name, maps in mapping.items():
+            if isinstance(maps, dict):
+                maps = [maps]
+            if any(
+                np.all(attval == var.attrs.get(att))  # attribute value is found in map
+                for ma in maps
+                for att, attval in ma.items()
+                if att in var.attrs  # variable has attribute
+            ):
+                # Must check that height is 10m or surface
+                for coordname, coordvar in var.coords.items():
+                    if coordvar.attrs.get('_CoordinateAxisType', None) == 'Height' and 10 in coordvar.values:
+                        logger.debug(f'Selecting GRIB variable {var_name} at 10m height and adding standard_name {standard_name}')
+                        ds[var_name] = var.sel({coordname: 10})
+                        ds[var_name].attrs['standard_name'] = standard_name
+    return ds
+
 def datetime_from_variable(var):
     import pandas as pd
     try:

opendrift/readers/reader_netCDF_CF_generic.py

@@ -24,7 +24,8 @@
 import pandas as pd
 import xarray as xr
 from opendrift.readers.basereader import BaseReader, StructuredReader
-from opendrift.readers import open_dataset_opendrift, datetime_from_variable
+from opendrift.readers import open_dataset_opendrift, datetime_from_variable, \
+                              add_standard_name_for_surface_grib_variables
 
 
 class Reader(StructuredReader, BaseReader):
@@ -115,6 +116,8 @@ def __init__(self, filename=None, zarr_storage_options=None, name=None, proj4=No
 
         self.Dataset = open_dataset_opendrift(source=filename, zarr_storage_options=zarr_storage_options)
 
+        self.Dataset = add_standard_name_for_surface_grib_variables(self.Dataset)
+
         if name is None:
             self.name = str(filename)
         else:
@@ -207,7 +210,12 @@ def __init__(self, filename=None, zarr_storage_options=None, name=None, proj4=No
             if standard_name == 'time' or axis == 'T' or var_name in ['time', 'vtime']:
                 # Read and store time coverage (of this particular file)
                 if len(var.dims)==1:
-                    self.dimensions['time'] = var.dims[0]
+                    if 'time' in self.dimensions:
+                        if not isinstance(self.dimensions['time'], list):
+                            self.dimensions['time'] = [self.dimensions['time']]
+                        self.dimensions['time'].append(var.dims[0])
+                    else:
+                        self.dimensions['time'] = var.dims[0]
 
                 self.times = datetime_from_variable(var)
                 if len(self.times) > 1:
@@ -378,6 +386,18 @@ def __init__(self, filename=None, zarr_storage_options=None, name=None, proj4=No
                     var = var.isel(ensemble_member=0).squeeze()
                     self.Dataset[va] = var
 
+        if isinstance(self.dimensions.get('time', None), list):
+            warning = f'Several time dimensions detected: {self.dimensions["time"]}'
+            variable_time_dimensions = []  # Check which time dimensions are actually used
+            for vn, va in self.variable_mapping.items():
+                variable_time_dimensions += [x for x in self.Dataset[va].coords if x in self.dimensions['time']]
+                variable_time_dimensions = list(set(variable_time_dimensions))
+            if len(variable_time_dimensions) == 1:
+                logger.debug(f'{warning}: Using single time dimension found in detected variables: {variable_time_dimensions}')
+                self.dimensions['time'] = variable_time_dimensions[0]
+            elif len(variable_time_dimensions) > 1:
+                raise ValueError(f'{warning}, and also several time dimensions used by detected variables: {variable_time_dimensions}')
+
         # Run constructor of parent Reader class
         super().__init__()
 

tests/models/test_models.py

@@ -70,7 +70,7 @@ def test_windblow(self):
         o.seed_elements(lon, lat, radius=5000, number=1000,
                         time=reader_arome.start_time)
         o.run(steps=24, time_step=3600)
-        self.assertAlmostEqual(o.elements.lon.max(), 15.864, 2)
+        self.assertAlmostEqual(o.elements.lon.max(), 15.881, 2)
 
     def test_shipdrift(self):
         """Sintef case study"""