Add CHELSA-ERA5 tiled climatologies.

hyoga/open/aggregator.py

@@ -9,6 +9,8 @@
 """
 
 import os.path
+import subprocess
+import tempfile
 import xarray as xr
 import hyoga.open.downloader
 
@@ -113,12 +115,87 @@ def pattern(self, *args):
         raise NotImplementedError("This should be implemented in subclasses.")
 
 
+class CERA5TiledAggregator(TiledAggregator):
+    """An aggregator to split CHELSA-ERA5 climatologies into tiles.
+
+    Call parameters
+    ---------------
+    variable : 'pr', 'tas', 'tasmax', 'tasmin'
+        The short name for the CHELSA-ERA5 variable aggregated among:
+        - daily mean precipitation ('pr', kg m-2 month-1),
+        - daily mean near-surface air temperature ('tas', K),
+        - daily maximum near surface air temperature ('tasmax', K),
+        - daily minimum near surface air temperature ('tasmin', K).
+    """
+
+    def inputs(self, *args):
+        """Return paths of input files, downloading as necessary."""
+        # FIXME consider moving tifs to e.g. .cache/hyoga/cera5/geotiff
+        variable, = args
+        downloader = hyoga.open.downloader.CacheDownloader()
+        basenames = [
+            f'CHELSA_{variable}_{month:02d}_1981-2010_V.2.1.tif'
+            for month in range(1, 13)]
+        paths = [downloader(
+            'https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/'
+            f'GLOBAL/climatologies/1981-2010/{variable}/{basename}',
+            f'chelsa/{basename}') for basename in basenames]
+        return paths
+
+    def pattern(self, *args):
+        variable, = args
+        xdg_cache = os.environ.get("XDG_CACHE_HOME", os.path.join(
+            os.path.expanduser('~'), '.cache'))
+        return os.path.join(
+            xdg_cache, 'hyoga', 'cera5', 'clim',
+            f'cera5.{variable}.mon.8110.avg.{{}}.nc')
+
+    def aggregate(self, inputs, output):
+        """Aggregate tiled `inputs` to files matching `output` pattern."""
+
+        # create directory if missing
+        os.makedirs(os.path.dirname(output[0]), exist_ok=True)
+
+        # open inputs as multi-file dataset
+        with xr.open_mfdataset(
+                inputs, combine='nested', concat_dim='month') as ds:
+
+            # for each tile
+            for tilepath, (lat, lon) in zip(output, self._get_all_coords()):
+
+                # check wether tile file exists
+                if os.path.isfile(tilepath):
+                    continue
+
+                # crop a 30x30 degree tile
+                # FIXME rename (x, y) to (lon, lat), band_data to (tas, pr)?
+                # FIXME flip decreasing latitudes to match CW5E5 data?
+                tile = ds.sel(x=slice(lon, lon+30), y=slice(lat+30, lat))
+                tile = tile.squeeze('band', drop=True)
+                tile.band_data.encoding.update(_FillValue=-9999)
+
+                # store as uncompressed netcdf
+                print(f"aggregating {tilepath} ...")
+                tile.to_netcdf(tilepath)
+
+                # nccopy compression beats xarray by far
+                print(f"compressing {tilepath} ...")
+                dirname, basename = os.path.split(tilepath)
+                with tempfile.NamedTemporaryFile(
+                        dir=dirname, prefix=basename+'.') as tmp:
+                    subprocess.run(['nccopy', '-sd6', tilepath, tmp.name])
+                    os.replace(tmp.name, tilepath)
+
+        # return multi-tile output pattern
+        return output
+
+
 class CW5E5TiledAggregator(TiledAggregator):
     """An aggregator to compute CHELSA-W5E5 climatologies from daily means.
 
     Call parameters
     ---------------
-    variable : 'tasmax', 'tas', 'tasmin', 'rsds', 'pr'
+    variable : 'pr', 'rsds', 'tas', 'tasmax', 'tasmin'
         The short name for the CHELSA-W5E5 variable aggregated among:
         - daily mean precipitation ('pr', kg m-2 s-1),
         - daily mean surface downwelling shortwave dadiation ('rsds', W m-2),

hyoga/open/reprojected.py

@@ -70,7 +70,14 @@ def _open_climatology(source='chelsa', variable='tas'):
             paths, combine='nested', concat_dim='time', decode_cf=True)
         da = ds.band_data.squeeze()
 
-    # CHELSA 1981-2010 global climatologies
+    # CHELSA-ERA5 1981-2010 global climatologies
+    elif source == 'cera5':
+        aggregator = hyoga.open.aggregator.CERA5TiledAggregator()
+        paths = aggregator(variable)
+        da = xr.open_mfdataset(paths, decode_coords='all').band_data
+        da.attrs.update(_FillValue=np.nan)  # so values show in ncview
+
+    # CHELSA-W5E5 1981-2010 global climatologies
     elif source == 'cw5e5':
         aggregator = hyoga.open.aggregator.CW5E5TiledAggregator()
         start, end = 1981, 2010  # FIXME allow custom aggregation period