From 2e5997c9d7a14111866c0f4c9eca71bc6ba92af6 Mon Sep 17 00:00:00 2001 From: ghiggi Date: Wed, 19 Jan 2022 15:39:56 +0100 Subject: [PATCH] Refactoring of upsampling utils + coverage --- pyresample/geometry.py | 193 +++++++++++++++++++++---------- pyresample/test/test_geometry.py | 57 +++++++++ 2 files changed, 192 insertions(+), 58 deletions(-) diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 45af21bd9..af76ea1c8 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -807,7 +807,7 @@ def upsample(self, x=1, y=1): swath_def.upsample(x=1, y=1) simply returns the current swath_def. """ # TODO: An alternative would be to use geotiepoints.geointerpolator.GeoInterpolator - # But I have some problem using it, see code snippet in the PR description. + # But I have some problem using it, see code snippet in a comment of the PR. import dask.array as da import numpy as np import pyproj @@ -824,70 +824,31 @@ def upsample(self, x=1, y=1): # Return SwathDefinition if nothing to upsample if x == 1 and y == 1: return self - - def _upsample_ranges_1D(x, factor=1): - ranges2D = np.linspace(x[:-1], x[1:], num=factor, endpoint=False, axis=1) - return np.concatenate((ranges2D.ravel(), [x[-1]])) - - def upsample_ranges_2D(x, factor=1, axis=0): - x = np.array(x) - if x.ndim not in [1, 2]: - raise ValueError("Expects 1D or 2D array.") - if not isinstance(axis, int): - raise TypeError("'axis' must be: 0 or 1 integer.") - if axis not in [0, 1]: - raise ValueError("Expects 'axis' 0 or 1") - if not isinstance(factor, int): - raise TypeError("'factor' must be an integer equal or larger than 1.") - if factor < 1: - raise ValueError("'factor' must be an integer equal or larger than 1.") - - if x.ndim == 1: - return _upsample_ranges_1D(x, factor=factor) - else: - l_ranges = [] - if axis == 1: - for i in range(x.shape[0]): - l_ranges.append(_upsample_ranges_1D(x[i, :], factor=factor)) - return np.vstack(l_ranges) - else: # axis = 0 - for i in range(x.shape[1]): - l_ranges.append(_upsample_ranges_1D(x[:, i], factor=factor)) - return np.vstack(l_ranges).transpose() - - def _upsample_corners(corners, x_factor=1, y_factor=1): - new_breaks_xx = upsample_ranges_2D(corners, factor=x_factor, axis=1) - new_corners = upsample_ranges_2D(new_breaks_xx, factor=y_factor, axis=0) - return new_corners + # --------------------------------------------------------------------. + # TODO: + # - Refactor for dask-compatibility + # - Should we make _infer_interval_breaks dask-compatible? def _get_corners_from_centroids(centroids): breaks_xx = _infer_interval_breaks(centroids, axis=1) corners = _infer_interval_breaks(breaks_xx, axis=0) return corners - def _get_centroids_from_corners(corners): - centroids = (corners[1:, 1:] + corners[:-1, :-1]) / 2 - return centroids - - # TODO: Decide if compute in memory or with dask - def upsample_centroids(centroid_x, centroid_y, centroid_z, x_factor=1, y_factor=1): - corners_x = _get_corners_from_centroids(centroid_x) - corners_y = _get_corners_from_centroids(centroid_y) - corners_z = _get_corners_from_centroids(centroid_z) - x_new_corners = _upsample_corners(corners_x, x_factor=x_factor, y_factor=y_factor) - y_new_corners = _upsample_corners(corners_y, x_factor=x_factor, y_factor=y_factor) - z_new_corners = _upsample_corners(corners_z, x_factor=x_factor, y_factor=y_factor) - x_new_centroids = _get_centroids_from_corners(x_new_corners) - y_new_centroids = _get_centroids_from_corners(y_new_corners) - z_new_centroids = _get_centroids_from_corners(z_new_corners) - return x_new_centroids, y_new_centroids, z_new_centroids - - def _upsample_centroid(centroid, x_factor=1, y_factor=1): + def _upsample_centroid(centroid, x=1, y=1): corners = _get_corners_from_centroids(centroid) - new_corners = _upsample_corners(corners, x_factor=x_factor, y_factor=y_factor) - new_centroids = _get_centroids_from_corners(new_corners) + # Retrieve corners of the the upsampled grid + new_corners = _linspace2D_between_values(corners, num_x=x - 1, num_y=y - 1) + # Get centroids from corners + new_centroids = (new_corners[:-1, :-1] + new_corners[1:, 1:]) / 2 return new_centroids + def upsample_centroids(centroid_x, centroid_y, centroid_z, x=1, y=1): + x_new_centroids = _upsample_centroid(centroid_x, x=x, y=y) + y_new_centroids = _upsample_centroid(centroid_y, x=x, y=y) + z_new_centroids = _upsample_centroid(centroid_z, x=x, y=y) + return x_new_centroids, y_new_centroids, z_new_centroids + + # --------------------------------------------------------------------. # Define geodetic and geocentric projection geocent = pyproj.Proj(proj='geocent') latlong = pyproj.Proj(proj='latlong') @@ -913,10 +874,9 @@ def _upsample_centroid(centroid, x_factor=1, y_factor=1): # x, # y) # res1 = xr.DataArray(res1, dims=['y', 'x', 'coord'], coords=src_lons.coords) - res = np.stack(upsample_centroids(res[:, :, 0].data, res[:, :, 1].data, - res[:, :, 2].data, x_factor=x, y_factor=y), axis=2) + res[:, :, 2].data, x=x, y=y), axis=2) new_centroids = xr.DataArray(da.from_array(res), dims=['y', 'x', 'xyz']) # Back-conversion to geographic CRS @@ -1267,6 +1227,123 @@ def _convert_2D_array(arr, to, dims=None): raise NotImplementedError +def _linspace1D_between_values(arr, num=0): + """Dask-friendly function linearly interpolating values between each 1D array values. + + This function does not perform extrapolation. + It expects a 1D array as input! + + Parameters + ---------- + arr : (np.ndarray, dask.array.Array) + Numpy or Dask Array to be linearly interpolated between values. + num : int, optional + The number of linearly spaced values to infer between array values. + The default is 0. + + Returns + ------- + arr : (np.ndarray, dask.array.Array) + Numpy or Dask Array with in-between linearly interpolated values. + + Example + ------- + + Function call: _linspace1D_between_values(arr, num=1) + Input array: + np.array([5.0, 7.0]) + + Output array: + np.array([5.0, 6.0, 7.0]) + """ + import xarray as xr + + # Check input validity + if arr.ndim != 1: + raise ValueError("'_linspace1D_between_values' expects a 1D array.") + num = int(num) + if num < 0: + raise ValueError("'x' and 'y' must be an integer equal or larger than 0.") + if num == 0: + return arr + # Define src and dst ties + dst_N = (arr.size - 1) * (num + 1) + 1 + src_ties = np.arange(dst_N, step=num + 1) + dst_ties = np.arange(dst_N) + # Interpolate + da = xr.DataArray( + data=arr, + dims=("x"), + coords={"x": src_ties} + ) + da_interp = da.interp(x=dst_ties, method="linear") + return da_interp.data + + +def _linspace2D_between_values(arr, num_x=0, num_y=0): + """Dask-friendly function linearly interpolating values between each 2D array values. + + This function does not perform extrapolation. + It expects a 2D array as input! + + Parameters + ---------- + arr : (np.ndarray, dask.array.Array) + Numpy or Dask Array to be linearly interpolated between values. + num_x : int, optional + The number of linearly spaced values to infer between array values (along x). + . The default is 0. + num_y : int, optional + The number of linearly spaced values to infer between array values (along y). + The default is 0. + + Returns + ------- + arr : (np.ndarray, dask.array.Array) + Numpy or Dask Array with in-between linearly interpolated values. + + Example + ------- + + Function call: _linspace1D_between_values(arr, num=1) + Input: + np.array([[5.0, 7.0], + [7.0, 9.0]]) + Output: + np.array([[5.0, 6.0, 7.0], + [6.0, 7.0, 8.0], + [7.0, 8.0, 9.0]]) + """ + import xarray as xr + + # Check input validity + if arr.ndim != 2: + raise ValueError("'_linspace2D_between_values' expects a 2D array.") + num_x = int(num_x) + num_y = int(num_y) + if num_x < 0 or num_y < 0: + raise ValueError("'x' and 'y' must be an integer equal or larger than 0.") + if num_x == 0 and num_y == 0: + return arr + # Define src and dst ties + shape = arr.shape + Nx_dst = (shape[1] - 1) * (num_x + 1) + 1 + Ny_dst = (shape[0] - 1) * (num_y + 1) + 1 + + src_ties_x = np.arange(Nx_dst, step=num_x + 1) + src_ties_y = np.arange(Ny_dst, step=num_y + 1) + dst_ties_x = np.arange(Nx_dst) + dst_ties_y = np.arange(Ny_dst) + # Interpolate + da = xr.DataArray( + data=arr, + dims=("y", "x"), + coords={"y": src_ties_y, "x": src_ties_x} + ) + da_interp = da.interp({"y": dst_ties_y, "x": dst_ties_x}, method="linear") + return da_interp.data + + def _get_extended_lonlats(lon_start, lat_start, lon_end, lat_end, npts, transpose=True): """Utils employed by SwathDefinition.extend. diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py index 27fb6115e..c842f4f00 100644 --- a/pyresample/test/test_geometry.py +++ b/pyresample/test/test_geometry.py @@ -2151,6 +2151,63 @@ def test_latslons_arr_conversion(self): self.assertRaises(TypeError, _convert_2D_array, [dict_format['Numpy']], 'numpy') self.assertRaises(ValueError, _convert_2D_array, dict_format['Numpy'], 'unvalid_format') + def test_linspace1D_between_values(self): + """Test linspace1D_between_values.""" + import dask.array as da + import numpy as np + + from pyresample.geometry import _linspace1D_between_values + arr_np = np.array([5.0, 7.0, 9.0]) + arr_dask = da.from_array(arr_np) + + res_np = _linspace1D_between_values(arr_np, num=1) + res_dask = _linspace1D_between_values(arr_dask, num=1) + + np.testing.assert_allclose(res_np, [5., 6., 7., 8., 9.]) + np.testing.assert_allclose(res_dask, [5., 6., 7., 8., 9.]) + assert isinstance(res_np, np.ndarray) + assert isinstance(res_dask, da.Array) + + # Test for no interpolation inbetween values + res = _linspace1D_between_values(arr_np, num=0) + np.testing.assert_allclose(res, arr_np) + + # Test for valid inputs + self.assertRaises(ValueError, _linspace1D_between_values, arr_np, -1) + self.assertRaises(ValueError, _linspace1D_between_values, np.zeros((2, 2)), 0) + + def test_linspace2D_between_values(self): + """Test linspace2D_between_values.""" + import dask.array as da + import numpy as np + + from pyresample.geometry import _linspace2D_between_values + arr_np = np.array([[5.0, 7.0], + [7.0, 9.0]]) + arr_dask = da.from_array(arr_np) + + res_np = _linspace2D_between_values(arr_np, num_x=1, num_y=3) + res_dask = _linspace2D_between_values(arr_dask, num_x=1, num_y=3) + + output_expected = np.array([[5., 6., 7.], + [5.5, 6.5, 7.5], + [6., 7., 8.], + [6.5, 7.5, 8.5], + [7., 8., 9.]]) + np.testing.assert_allclose(res_np, output_expected) + np.testing.assert_allclose(res_dask, output_expected) + assert isinstance(res_np, np.ndarray) + assert isinstance(res_dask, da.Array) + + # Test for no interpolation inbetween values + res = _linspace2D_between_values(arr_np, num_x=0, num_y=0) + np.testing.assert_allclose(res, arr_np) + + # Test for valid inputs + self.assertRaises(ValueError, _linspace2D_between_values, arr_np, -1, 0) + self.assertRaises(ValueError, _linspace2D_between_values, arr_np, 0, -1) + self.assertRaises(ValueError, _linspace2D_between_values, arr_np[0, :], 0, 0) + def test_get_extended_lonlats(self): import numpy as np