Add extra steps to diagnostic to make output of hydrology/recipe_lisflood.yml usable by the LISFLOOD model (ESMValGroup#1737)

* draft of extra steps to make output forcing usable by lisflood

* replace draft cube manipulations with iris in diagnostic and process data over global domain in recipe

* clean up and turn crop region

* loop over coordinates to remove unused coordinates and the bounds of used coordinates

* reorder coordinates

* change to xarray as sorting variables and coordinates is not possible in iris

* fix output file encoding

* fix wrong variable name

* fix date fetching

* Apply suggestions from code review

* Automatic formatting

Co-authored-by: Peter Kalverla <peter.kalverla@gmx.com>

Author: JaroCamphuijsen

esmvaltool/diag_scripts/hydrology/lisflood.py

@@ -3,6 +3,8 @@
 from pathlib import Path
 
 import iris
+import numpy as np
+import xarray as xr
 from iris.analysis.maths import exp as iris_exp
 
 from esmvaltool.diag_scripts.shared import (ProvenanceLogger,
@@ -20,6 +22,7 @@ def get_provenance_record(ancestor_files):
         'authors': [
             'verhoeven_stefan',
             'kalverla_peter',
+            'camphuijsen_jaro',
         ],
         'projects': [
             'ewatercycle',
@@ -76,16 +79,15 @@ def compute_vapour_pressure(tdps):
     esat.long_name = 'Daily Actual Water Vapour Pressure'
     esat.standard_name = 'water_vapor_pressure'
     esat.units = 'hPa'
-    esat.attributes['comment'] = ''.join((
-        'Actual water vapour pressure of air near the surface calculated',
-        ' from tdps using Tetens formula'
-    ))
+    esat.attributes['comment'] = ''.join(
+        ('Actual water vapour pressure of air near the surface calculated',
+         ' from tdps using Tetens formula'))
     return esat
 
 
 def compute_windspeed(uas, vas):
     """Compute absolute wind speed from horizontal components."""
-    sfc_wind = (uas ** 2 + vas ** 2) ** .5
+    sfc_wind = (uas**2 + vas**2)**.5
     sfc_wind.var_name = 'sfcWind'
     sfc_wind.long_name = 'Daily-Mean Near-Surface Wind Speed'
     sfc_wind.standard_name = 'wind_speed'
@@ -94,11 +96,10 @@ def compute_windspeed(uas, vas):
     return sfc_wind
 
 
-def save(cube, var_name, dataset, cfg):
+def save(xrds, var_name, dataset, cfg):
     """Save processed cube to a lisflood-compatible file."""
-    time_coord = cube.coord('time')
-    start_year = time_coord.cell(0).point.year
-    end_year = time_coord.cell(-1).point.year
+    start_year = int(xrds.time[0].dt.year)
+    end_year = int(xrds.time[-1].dt.year)
     basename = '_'.join([
         'lisflood',
         dataset,
@@ -108,7 +109,7 @@ def save(cube, var_name, dataset, cfg):
         str(end_year),
     ])
     output_file = get_diagnostic_filename(basename, cfg)
-    iris.save(cube, output_file, fill_value=1.e20)
+    xrds.to_netcdf(output_file, encoding={var_name: {'_FillValue': 1.e20}})
     return output_file
 
 
@@ -138,19 +139,23 @@ def main(cfg):
         cubes['pr'].units = 'mm d-1'
 
         for var_name, cube in cubes.items():
-            cube.remove_coord('shape_id')
             # Western emisphere longitudes should be negative
             points = cube.coord('longitude').points
             cube.coord('longitude').points = (points + 180) % 360 - 180
             # latitudes decreasing
             cube = cube[:, ::-1, ...]
 
-            output_file = save(cube, var_name, dataset, cfg)
+            # convert to xarray dataset (xrds)
+            # remove coordinate bounds drop extra coordinates and reorder
+            xrds = xr.DataArray.from_iris(cube).to_dataset()
+            ordered_coords = ['lon', 'lat', 'time']
+            extra_coords = np.setdiff1d(xrds.coords, ordered_coords)
+            xrds = xrds.drop(extra_coords)[ordered_coords + [var_name]]
+
+            output_file = save(xrds, var_name, dataset, cfg)
 
             # Store provenance
-            provenance_record = get_provenance_record(
-                ancestors[var_name]
-            )
+            provenance_record = get_provenance_record(ancestors[var_name])
             with ProvenanceLogger(cfg) as provenance_logger:
                 provenance_logger.log(output_file, provenance_record)