Fix bilinear resampling to areas with invalid coordinates

pyresample/bilinear/_base.py

@@ -143,13 +143,15 @@ def _create_empty_bil_info(self):
         self.mask_slices = self._index_array >= self._source_geo_def.size
 
     def _get_valid_output_indices(self):
-        return ((self._target_lons >= -180) & (self._target_lons <= 180) &
-                (self._target_lats <= 90) & (self._target_lats >= -90))
+        self._valid_output_indices = np.ravel(
+            (self._target_lons >= -180) & (self._target_lons <= 180) &
+            (self._target_lats <= 90) & (self._target_lats >= -90))
 
     def _get_index_array(self):
+        self._get_valid_output_indices()
         index_array = _query_no_distance(
             self._target_lons, self._target_lats,
-            self._get_valid_output_indices(), self._resample_kdtree,
+            self._valid_output_indices, self._resample_kdtree,
             self._neighbours, self._epsilon,
             self._radius_of_influence)
         self._index_array = self._reduce_index_array(index_array)
@@ -170,7 +172,10 @@ def _get_fractional_distances(self):
             corner_points, out_x, out_y)
 
     def _get_output_xy(self):
-        return _get_output_xy(self._target_geo_def)
+        out_x, out_y = _get_output_xy(self._target_geo_def)
+        out_x = out_x[self._valid_output_indices]
+        out_y = out_y[self._valid_output_indices]
+        return out_x, out_y
 
     def _get_input_xy(self):
         return _get_input_xy(self._source_geo_def,
@@ -637,9 +642,8 @@ def _resample(corner_points, fractional_distances):
 def _query_no_distance(target_lons, target_lats,
                        valid_output_index, kdtree, neighbours, epsilon, radius):
     """Query the kdtree. No distances are returned."""
-    voir = np.ravel(valid_output_index)
-    target_lons_valid = np.ravel(target_lons)[voir]
-    target_lats_valid = np.ravel(target_lats)[voir]
+    target_lons_valid = np.ravel(target_lons)[valid_output_index]
+    target_lats_valid = np.ravel(target_lats)[valid_output_index]
 
     _, index_array = kdtree.query(
         lonlat2xyz(target_lons_valid, target_lats_valid),

pyresample/bilinear/xarr.py

@@ -88,7 +88,10 @@ def _get_input_xy(self):
                              self._valid_input_index, self._index_array)
 
     def _get_output_xy(self):
-        return _get_output_xy(self._target_geo_def)
+        out_x, out_y = _get_output_xy(self._target_geo_def)
+        out_x = out_x[self._valid_output_indices]
+        out_y = out_y[self._valid_output_indices]
+        return out_x, out_y
 
     def _limit_output_values_to_input(self, data, res, fill_value):
         epsilon = 1e-6
@@ -102,6 +105,19 @@ def _limit_output_values_to_input(self, data, res, fill_value):
         return da.where(np.isnan(res), fill_value, res)
 
     def _reshape_to_target_area(self, res, ndim):
+        if ndim == 3:
+            dim_multiplier = res.shape[0]
+        else:
+            dim_multiplier = 1
+            res = da.reshape(res, (1, res.size))
+        if res.size != dim_multiplier * self._target_geo_def.size:
+            out = []
+            for i in range(dim_multiplier):
+                tmp = da.full(self._target_geo_def.size, np.nan)
+                tmp[self._valid_output_indices] = res[i, :]
+                out.append(tmp)
+            res = da.stack(out)
+
         shp = self._target_geo_def.shape
         if ndim == 3:
             res = da.reshape(res, (res.shape[0], shp[0], shp[1]))

pyresample/test/test_bilinear.py

@@ -1032,7 +1032,8 @@ def test_query_no_distance(self):
         kdtree = mock.MagicMock()
         kdtree.query.return_value = (1, 2)
         lons, lats = self.target_def.get_lonlats()
-        voi = (lons >= -180) & (lons <= 180) & (lats <= 90) & (lats >= -90)
+        voi = np.ravel(
+            (lons >= -180) & (lons <= 180) & (lats <= 90) & (lats >= -90))
         res = _query_no_distance(lons, lats, voi, kdtree, self._neighbours,
                                  0., self.radius)
         # Only the second value from the query is returned
@@ -1159,3 +1160,47 @@ def test_check_fill_value():
 
     # float fill value + float dtype -> no change
     assert _check_fill_value(3.3, np.float32)
+
+
+def test_target_has_invalid_coordinates():
+    """Test bilinear resampling to area that has invalid coordinates.
+
+    The area used here is in geos projection that has space pixels in the corners.
+    """
+    import dask.array as da
+    import xarray as xr
+
+    from pyresample.bilinear import XArrayBilinearResampler
+
+    # NumpyBilinearResampler
+    from pyresample.geometry import AreaDefinition, GridDefinition
+
+    geos_def = AreaDefinition('geos',
+                              'GEO area with space in corners',
+                              'geos',
+                              {'proj': 'geos',
+                               'lon_0': '0.0',
+                               'a': '6378169.0',
+                               'b': '6356583.8',
+                               'h': '35785831.0'},
+                              640, 640,
+                              [-5432229.931711678,
+                               -5429229.528545862,
+                               5429229.528545862,
+                               5432229.931711678])
+    lats = np.linspace(-89, 89, 179)
+    lons = np.linspace(-179, 179, 359)
+    lats = np.repeat(lats[:, None], 359, axis=1)
+    lons = np.repeat(lons[None, :], 179, axis=0)
+    grid_def = GridDefinition(lons=lons, lats=lats)
+
+    data_xr = xr.DataArray(da.random.uniform(0, 1, lons.shape), dims=["y", "x"])
+
+    resampler = XArrayBilinearResampler(grid_def,
+                                        geos_def,
+                                        500e3,
+                                        reduce_data=False)
+    res = resampler.resample(data_xr)
+    res = res.compute()
+    assert not np.all(np.isnan(res))
+    assert np.any(np.isnan(res))