From 1f086ff737da8a2faef09421b43f2c4bba56b5ba Mon Sep 17 00:00:00 2001 From: Norman Fomferra Date: Mon, 30 Apr 2018 19:01:48 +0200 Subject: [PATCH] Support datasets with 0,360 degree longitude ranges closes #620 --- CHANGES.md | 1 + cate/core/cdm.py | 9 +++++-- cate/ops/__init__.py | 2 +- cate/ops/normalize.py | 43 +++++++++++++++++++++++++++++++++- cate/util/im/geoextent.py | 4 +++- test/ops/test_normalize.py | 26 +++++++++++++++++++- test/util/im/test_geoextent.py | 5 ++++ 7 files changed, 84 insertions(+), 6 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index eddae5709..412388509 100644 --- a/CHANGES.md +++ b/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) diff --git a/cate/core/cdm.py b/cate/core/cdm.py index 2ebb8a3e6..306b735e8 100644 --- a/cate/core/cdm.py +++ b/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: diff --git a/cate/ops/__init__.py b/cate/ops/__init__.py index 305f60116..761602381 100644 --- a/cate/ops/__init__.py +++ b/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) diff --git a/cate/ops/normalize.py b/cate/ops/normalize.py index 798b4b4be..4e848a518 100644 --- a/cate/ops/normalize.py +++ b/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 diff --git a/cate/util/im/geoextent.py b/cate/util/im/geoextent.py index a4c0bc9ab..1af5fd998 100644 --- a/cate/util/im/geoextent.py +++ b/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: diff --git a/test/ops/test_normalize.py b/test/ops/test_normalize.py index f35336e59..be6148cd5 100644 --- a/test/ops/test_normalize.py +++ b/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]) diff --git a/test/util/im/test_geoextent.py b/test/util/im/test_geoextent.py index 6e6a3bd07..6acf6590d 100644 --- a/test/util/im/test_geoextent.py +++ b/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