Merge pull request #249 from sfinkens/bucket-lowres

Fix bucket assignment

pyresample/bucket/__init__.py

@@ -90,24 +90,17 @@ 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
         ----------
         lons : Numpy or Dask array
             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) &

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."""