Support datasets with 0,360 degree longitude ranges

closes #620

CHANGES.md

@@ -1,5 +1,6 @@
 ## Version 2.0.0.dev10 (in development)
 
+* Support datasets with 0,360 degree longitude ranges [#620](https://github.com/CCI-Tools/cate/issues/620)
 * Temporal aggregation operation can now aggregate to pre-defined seasons, as well as custom resolutions [#472](https://github.com/CCI-Tools/cate/issues/472)
 * We now use "MB" units instead of "MiB" (part of [#325](https://github.com/CCI-Tools/cate/issues/325))
 * Fixed a bug with animation generation [#585](https://github.com/CCI-Tools/cate/issues/585)

cate/core/cdm.py

@@ -107,7 +107,7 @@
 Components
 ==========
 """
-
+import warnings
 from collections import OrderedDict
 from typing import List, Optional, Union
 
@@ -355,7 +355,12 @@ def get_tiling_scheme(var: xr.DataArray) -> Optional[TilingScheme]:
     width, height = var.shape[-1], var.shape[-2]
     lats = var.coords[lat_dim_name]
     lons = var.coords[lon_dim_name]
-    geo_extent = GeoExtent.from_coord_arrays(lons, lats)
+    try:
+        geo_extent = GeoExtent.from_coord_arrays(lons, lats)
+    except ValueError as e:
+        warnings.warn(f'failed to derive geo-extent for tiling scheme: {e}')
+        # Create a default geo-extent which is probably wrong, but at least we see something
+        geo_extent = GeoExtent()
     try:
         return TilingScheme.create(width, height, 360, 360, geo_extent)
     except ValueError:

cate/ops/__init__.py

@@ -51,7 +51,7 @@ def cate_init():
 from .select import select_var
 from .coregistration import coregister
 from .correlation import pearson_correlation_scalar, pearson_correlation
-from .normalize import normalize, adjust_temporal_attrs, adjust_spatial_attrs
+from .normalize import normalize, adjust_temporal_attrs, adjust_spatial_attrs, fix_lon_360
 from .io import (open_dataset, save_dataset, read_object, write_object,
                  read_text, write_text, read_json, write_json, read_csv,
                  read_geo_data_frame, read_netcdf, write_netcdf3, write_netcdf4)

cate/ops/normalize.py

@@ -30,9 +30,11 @@
 """
 
 import xarray as xr
+import numpy as np
 
 from cate.core.op import op, op_return
 from cate.core.opimpl import normalize_impl, adjust_temporal_attrs_impl, adjust_spatial_attrs_impl
+from cate.core.types import ValidationError
 
 
 @op(tags=['utility'], version='1.0')
@@ -60,7 +62,7 @@ def normalize(ds: xr.Dataset) -> xr.Dataset:
 
 
 @op(tags=['utility'], version='1.0')
-def adjust_spatial_attrs(ds: xr.Dataset, allow_point: bool=False) -> xr.Dataset:
+def adjust_spatial_attrs(ds: xr.Dataset, allow_point: bool = False) -> xr.Dataset:
     """
     Adjust the global spatial attributes of the dataset by doing some
     introspection of the dataset and adjusting the appropriate attributes
@@ -101,3 +103,42 @@ def adjust_temporal_attrs(ds: xr.Dataset) -> xr.Dataset:
     :return: Adjusted dataset
     """
     return adjust_temporal_attrs_impl(ds)
+
+
+@op(tags=['utility'], version='1.0')
+def fix_lon_360(ds: xr.Dataset) -> xr.Dataset:
+    """
+    Fix the longitude of the given dataset ``ds`` so that it ranges from -180 to +180 degrees.
+
+    :param ds: The dataset whose longitudes are given in the range 0 to 360.
+    :return: The fixed dataset.
+    """
+    if 'lon' not in ds.coords:
+        raise ValidationError('missing coordinate variable "lon"')
+    if 'lon' not in ds.sizes:
+        raise ValidationError('missing dimension "lon"')
+    if len(ds.lon.shape) != 1:
+        raise ValidationError('coordinate variable "lon" must be 1-dimensional')
+    if len(ds.lon) < 2:
+        raise ValidationError('coordinate variable "lon" must have more than one element')
+
+    new_ds = ds.copy()
+    lon_size = ds.sizes['lon']
+    lon_size_05 = lon_size // 2
+
+    for var_name in new_ds.variables:
+        if var_name != 'lon':
+            var = new_ds.variables[var_name]
+            if len(var.dims) >= 1 and var.dims[-1] == 'lon':
+                temp = var.values[..., : lon_size_05]
+                var.values[..., : lon_size_05] = var.values[..., lon_size_05:]
+                var.values[..., lon_size_05:] = temp
+
+    delta_lon = new_ds['lon'][1] - new_ds['lon'][0]
+
+    new_ds['lon'] = xr.DataArray(np.linspace(-180. + 0.5 * delta_lon, +180. - 0.5 * delta_lon, lon_size),
+                                 dims=ds['lon'].dims,
+                                 attrs=ds['lon'].attrs)
+
+    new_ds['lon'].attrs['units'] = 'degrees east'
+    return new_ds

cate/util/im/geoextent.py

@@ -138,6 +138,7 @@ def from_coord_arrays(cls, x: np.ndarray, y: np.ndarray, eps: float = EPS) -> 'G
             y = y[(0,) * (y.ndim - 2) + (..., 0)]
 
         dx = None
+        x_norm = False
         if x.size > 1:
             dx = np.gradient(x)
             if (dx.max() - dx.min()) >= eps:
@@ -146,6 +147,7 @@ def from_coord_arrays(cls, x: np.ndarray, y: np.ndarray, eps: float = EPS) -> 'G
                 if x[0] > x[-1]:
                     # normalize x
                     x = np.where(x < 0., 360. + x, x)
+                    x_norm = True
                     # and test once more
                     dx = np.gradient(x)
                     fail = (dx.max() - dx.min()) >= eps
@@ -174,7 +176,7 @@ def from_coord_arrays(cls, x: np.ndarray, y: np.ndarray, eps: float = EPS) -> 'G
         y1 = y[0] - 0.5 * dy
         y2 = y[-1] + 0.5 * dy
 
-        if x2 > 180.0:
+        if x_norm and x2 > 180.0:
             x2 -= 360.0
 
         if y1 < y2:

test/ops/test_normalize.py

@@ -9,9 +9,10 @@
 import numpy as np
 import xarray as xr
 from jdcal import gcal2jd
+from numpy.testing import assert_array_almost_equal
 
 from cate.core.op import OP_REGISTRY
-from cate.ops.normalize import normalize, adjust_spatial_attrs, adjust_temporal_attrs
+from cate.ops.normalize import normalize, adjust_spatial_attrs, adjust_temporal_attrs, fix_lon_360
 from cate.util.misc import object_to_qualified_name
 
 
@@ -606,3 +607,26 @@ def test_only_time_dim_generated(self):
         self.assertEqual(new_ds.first.shape, (1, 90, 180))
         self.assertEqual(new_ds.second.shape, (1, 90, 180))
         self.assertEqual(new_ds.coords['time'][0], xr.DataArray(pd.to_datetime('2012-01-01')))
+
+
+class TestFix360Lon(TestCase):
+    def test_fix_360_lon(self):
+        lon_size = 360
+        lat_size = 130
+        time_size = 12
+        ds = xr.Dataset({
+            'first': (['time', 'lat', 'lon'], np.random.random_sample([time_size, lat_size, lon_size])),
+            'second': (['time', 'lat', 'lon'], np.random.random_sample([time_size, lat_size, lon_size])),
+            'lon': np.linspace(1., 360., lon_size),
+            'lat': np.linspace(-65, 65, lat_size),
+            'time': [datetime(2000, x, 1) for x in range(1, time_size + 1)]})
+
+        new_ds = fix_lon_360(ds)
+        self.assertIsNot(ds, new_ds)
+        self.assertEqual(ds.dims, new_ds.dims)
+        self.assertEqual(ds.sizes, new_ds.sizes)
+        assert_array_almost_equal(new_ds.lon, np.linspace(-179.5, 179.5, 360))
+        assert_array_almost_equal(new_ds.first[..., :180], ds.first[..., 180:])
+        assert_array_almost_equal(new_ds.first[..., 180:], ds.first[..., :180])
+        assert_array_almost_equal(new_ds.second[..., :180], ds.second[..., 180:])
+        assert_array_almost_equal(new_ds.second[..., 180:], ds.second[..., :180])

test/util/im/test_geoextent.py

@@ -74,6 +74,11 @@ def test_from_coord_arrays(self):
         np.testing.assert_almost_equal(np.array(rect.coords), np.array((-180.0, -90.0, 180.0, 90.0)))
         self.assertEqual(rect.inv_y, False)
 
+        # Example from https://research.csiro.au/slrwavescoast/sea-level/measurements-and-data/sea-level-data/
+        with self.assertRaises(ValueError) as cm:
+            GeoExtent.from_coord_arrays(np.arange(1., 361., 1.), np.arange(-65., 66., 1.))
+        self.assertEqual(str(cm.exception), 'east out of bounds: 360.5')
+
     def test_from_coord_arrays_with_eps(self):
         eps = 1e-4
         eps025 = 0.25 * eps