diff --git a/pyresample/bucket/__init__.py b/pyresample/bucket/__init__.py index 9e87678fc..9b9768610 100644 --- a/pyresample/bucket/__init__.py +++ b/pyresample/bucket/__init__.py @@ -90,8 +90,8 @@ def __init__(self, target_area, source_lons, source_lats): self._get_indices() self.counts = None - def _get_proj_coordinates(self, lons, lats, x_res, y_res): - """Calculate projection coordinates and round to resolution unit. + def _get_proj_coordinates(self, lons, lats): + """Calculate projection coordinates. Parameters ---------- @@ -99,15 +99,8 @@ def _get_proj_coordinates(self, lons, lats, x_res, y_res): Longitude coordinates lats : Numpy or Dask array Latitude coordinates - x_res : float - Resolution of the output in X direction - y_res : float - Resolution of the output in Y direction """ proj_x, proj_y = self.prj(lons, lats) - proj_x = round_to_resolution(proj_x, x_res) - proj_y = round_to_resolution(proj_y, y_res) - return np.stack((proj_x, proj_y)) def _get_indices(self): @@ -121,30 +114,20 @@ def _get_indices(self): Y indices of the target grid where the data are put """ LOG.info("Determine bucket resampling indices") - adef = self.target_area + # Transform source lons/lats to target projection coordinates x/y lons = self.source_lons.ravel() lats = self.source_lats.ravel() - - # Create output grid coordinates in projection units - x_res = (adef.area_extent[2] - adef.area_extent[0]) / adef.width - y_res = (adef.area_extent[3] - adef.area_extent[1]) / adef.height - x_vect = da.arange(adef.area_extent[0] + x_res / 2., - adef.area_extent[2] - x_res / 2., x_res) - # Orient so that 0-meridian is pointing down - y_vect = da.arange(adef.area_extent[3] - y_res / 2., - adef.area_extent[1] + y_res / 2., - -y_res) - - result = da.map_blocks(self._get_proj_coordinates, lons, - lats, x_res, y_res, + result = da.map_blocks(self._get_proj_coordinates, lons, lats, new_axis=0, chunks=(2,) + lons.chunks) proj_x = result[0, :] proj_y = result[1, :] - # Calculate array indices - x_idxs = ((proj_x - np.min(x_vect)) / x_res).astype(np.int) - y_idxs = ((np.max(y_vect) - proj_y) / y_res).astype(np.int) + # Calculate array indices. Orient so that 0-meridian is pointing down. + adef = self.target_area + x_res, y_res = adef.resolution + x_idxs = da.floor((proj_x - adef.area_extent[0]) / x_res).astype(np.int) + y_idxs = da.floor((adef.area_extent[3] - proj_y) / y_res).astype(np.int) # Get valid index locations mask = ((x_idxs >= 0) & (x_idxs < adef.width) & diff --git a/pyresample/test/test_bucket.py b/pyresample/test/test_bucket.py index bf1583944..c936bbfc3 100644 --- a/pyresample/test/test_bucket.py +++ b/pyresample/test/test_bucket.py @@ -5,6 +5,7 @@ import xarray as xr from unittest.mock import MagicMock, patch +from pyresample import create_area_def from pyresample.geometry import AreaDefinition from pyresample import bucket from pyresample.test.utils import CustomScheduler @@ -19,7 +20,6 @@ class Test(unittest.TestCase): 'lon_0': '0.0', 'proj': 'stere'}, 2560, 2048, (-3780000.0, -7644000.0, 3900000.0, -1500000.0)) - chunks = 2 lons = da.from_array(np.array([[25., 25.], [25., 25.]]), chunks=chunks) @@ -70,14 +70,12 @@ def test_get_proj_coordinates(self): prj.return_value = ([3.1, 3.1, 3.1], [4.8, 4.8, 4.8]) lons = [1., 1., 1.] lats = [2., 2., 2.] - x_res, y_res = 0.5, 0.5 self.resampler.prj = prj - result = self.resampler._get_proj_coordinates(lons, lats, x_res, y_res) + result = self.resampler._get_proj_coordinates(lons, lats) prj.assert_called_once_with(lons, lats) self.assertTrue(isinstance(result, np.ndarray)) - self.assertEqual(result.shape, (2, 3)) - self.assertTrue(np.all(result == np.array([[3., 3., 3.], - [5., 5., 5.]]))) + np.testing.assert_equal(result, np.array([[3.1, 3.1, 3.1], + [4.8, 4.8, 4.8]])) def test_get_bucket_indices(self): """Test calculation of array indices.""" @@ -86,8 +84,28 @@ def test_get_bucket_indices(self): self.resampler._get_indices() x_idxs, y_idxs = da.compute(self.resampler.x_idxs, self.resampler.y_idxs) - self.assertTrue(np.all(x_idxs == np.array([1709, 1709, 1706, 1705]))) - self.assertTrue(np.all(y_idxs == np.array([465, 465, 458, 455]))) + np.testing.assert_equal(x_idxs, np.array([1710, 1710, 1707, 1705])) + np.testing.assert_equal(y_idxs, np.array([465, 465, 459, 455])) + + # Additional small test case + adef = create_area_def( + area_id='test', + projection={'proj': 'latlong'}, + width=2, height=2, + center=(0, 0), + resolution=10) + lons = da.from_array( + np.array([-10.0, -9.9, -0.1, 0, 0.1, 9.9, 10.0, -10.1, 0]), + chunks=2) + lats = da.from_array( + np.array([-10.0, -9.9, -0.1, 0, 0.1, 9.9, 10.0, 0, 10.1]), + chunks=2) + resampler = bucket.BucketResampler(source_lats=lats, + source_lons=lons, + target_area=adef) + resampler._get_indices() + np.testing.assert_equal(resampler.x_idxs, np.array([-1, 0, 0, 1, 1, 1, -1, -1, -1])) + np.testing.assert_equal(resampler.y_idxs, np.array([-1, 1, 1, 1, 0, 0, -1, -1, -1])) def test_get_sum(self): """Test drop-in-a-bucket sum."""