Cannot resample with `bilinear` from lat/lon grid onto MSG full disk

[simonrp84]

Firstly, I'm not sure if this is a pyresample issue or a satpy issue, but I think the former is more likely so I'll put this here.

I'm trying to resample some gridded data onto the SEVIRI disk and am coming up against an error:

from pyresample.geometry import GridDefinition
from satpy.resample import get_area_def
from satpy import resample
import xarray as xr
import numpy as np

# Target area
targ_area = get_area_def('msg_seviri_fes_3km')

# Set up lats and lons
lats = np.linspace(-89, 89, 357)
lons = np.linspace(-179.75, 179.75, 720)
lats = np.repeat(lats[:, None], 720, axis=1)
lons = np.repeat(lons[None, :], 357, axis=0)

# Create grid definition
grid_def = GridDefinition(lons=lons, lats=lats)

# Create fake data
data_xr = xr.DataArray(np.random.uniform(0, 1, lons.shape), dims=["y", "x"])

# Do resampling
res = resample.resample(grid_def,
                        data_xr,
                        targ_area,
                        resampler='bilinear',
                        reduce_data=False,
                        radius_of_influence=50000)

Which results in the following error:

---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
File c:\users\simon\pycharmprojects\satpy\satpy\resample.py:873, in BilinearResampler.precompute(***failed resolving arguments***)
    872 try:
--> 873     self.load_bil_info(cache_dir, **kwargs)
    874     LOG.debug("Loaded bilinear parameters")

File c:\users\simon\pycharmprojects\satpy\satpy\resample.py:894, in BilinearResampler.load_bil_info(self, cache_dir, **kwargs)
    893 else:
--> 894     raise IOError

OSError: 

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Input In [23], in <module>
----> 1 res = resample.resample(grid_def,
      2                         data_xr,
      3                         targ_area,
      4                         resampler='bilinear',
      5                         reduce_data=False,
      6                         radius_of_influence=50000)

File c:\users\simon\pycharmprojects\satpy\satpy\resample.py:1358, in resample(source_area, data, destination_area, resampler, **kwargs)
   1356     res = [resampler_instance.resample(ds, **kwargs) for ds in data]
   1357 else:
-> 1358     res = resampler_instance.resample(data, **kwargs)
   1360 return res

File c:\users\simon\pycharmprojects\satpy\satpy\resample.py:438, in BaseResampler.resample(self, data, cache_dir, mask_area, **kwargs)
    435         kwargs['mask'] = data.isnull()
    436     kwargs['mask'] = kwargs['mask'].all(dim=flat_dims)
--> 438 cache_id = self.precompute(cache_dir=cache_dir, **kwargs)
    439 return self.compute(data, cache_id=cache_id, **kwargs)

File c:\users\simon\pycharmprojects\satpy\satpy\resample.py:877, in BilinearResampler.precompute(***failed resolving arguments***)
    875 except IOError:
    876     LOG.debug("Computing bilinear parameters")
--> 877     self.resampler.get_bil_info()
    878     LOG.debug("Saving bilinear parameters.")
    879     self.save_bil_info(cache_dir, **kwargs)

File ~\miniconda3\lib\site-packages\pyresample\bilinear\_base.py:118, in BilinearBase.get_bil_info(self, kdtree_class, nprocs)
    115 self._get_index_array()
    117 # Calculate vertical and horizontal fractional distances t and s
--> 118 self._get_fractional_distances()
    119 self._get_target_proj_vectors()
    120 self._get_slices()

File ~\miniconda3\lib\site-packages\pyresample\bilinear\_base.py:166, in BilinearBase._get_fractional_distances(self)
    163 out_x, out_y = self._get_output_xy()
    164 # Get the four closest corner points around each output location
    165 corner_points, self._index_array = \
--> 166     _get_four_closest_corners(*self._get_input_xy(),
    167                               out_x, out_y,
    168                               self._neighbours, self._index_array)
    169 self.bilinear_t, self.bilinear_s = _get_fractional_distances(
    170     corner_points, out_x, out_y)

File ~\miniconda3\lib\site-packages\pyresample\bilinear\_base.py:319, in _get_four_closest_corners(in_x, in_y, out_x, out_y, neighbours, index_array)
    310 """Get bounding corners.
    311 
    312 Get four closest locations from (in_x, in_y) so that they form a
   (...)
    315 
    316 """
    317 # Find four closest pixels around the target location
--> 319 stride, valid_corners = _get_stride_and_valid_corner_indices(
    320     out_x, out_y, in_x, in_y, neighbours)
    321 res = []
    322 indices = []

File ~\miniconda3\lib\site-packages\pyresample\bilinear\_base.py:521, in _get_stride_and_valid_corner_indices(out_x, out_y, in_x, in_y, neighbours)
    518 out_y_tile = _tile_output_coordinate_vector(out_y, neighbours)
    520 # Get differences in both directions
--> 521 x_diff = out_x_tile - in_x
    522 y_diff = out_y_tile - in_y
    524 return (np.arange(x_diff.shape[0]), (
    525     (x_diff > 0) & (y_diff < 0),  # Upper left corners
    526     (x_diff < 0) & (y_diff < 0),  # Upper right corners
    527     (x_diff > 0) & (y_diff > 0),  # Lower left corners
    528     (x_diff < 0) & (y_diff > 0))  # Lower right corners
    529 )

ValueError: operands could not be broadcast together with shapes (13778944,32) (10280821,32)

I've tried various combinations of radius_of_influence and reduce_data and always get the same error. Switching to nearest or gradient_search works, but I would prefer to use bilinear in this case as the results should be better quality for what I need. Any ideas what the problem is?

(edit) I also tried this with SwathDefinition instead of GridDefinition and I get the same error.

[pnuu]

The first thought was that it must be the space pixels, but no, the size 13778944 matches the size of the FD data (3712x3712). I'll have another look in the morning, getting too tired to figure out what the resampler is doing.

[pnuu]

Debugging was fun. The trace points to pyresample.bilinear._base, but some of the calls actually happen in .xarr which overloads some of the methods. The index_array used to expand (slice) the input X/Y coordinates to the output grid is too small to get full shape. I started to think that the resampler ever works only if the output grid is larger in pixels than the input data. But testing that by multiplying the input grid dimensions by 10 I got the exact same failure with the same mismatching shapes.

Need to do other things now, but will get back to this.

[simonrp84]

Thanks for looking at this Panu!

[pnuu]

Hah, it is the space pixels where both lons and lats are Inf. Via debugger:

ipdb> self._get_valid_output_indices().sum()
10280821
ipdb> self._target_lons.size - np.isinf(self._target_lons).sum()
10280821

If I mask these space coordinates to np.nan, the result is the same, so I think KDTree does something similar for nan coordinates. To further validate this hypothesis, I cropped the target area not to contain any space pixels by targ_area = get_area_def('msg_seviri_fes_3km')[1000:2000, 1000:2000] and that works.

So, how to fix this? Not sure. The tiling of output X and Y coordinates done in pyresample.bilinear._base._get_stride_and_valid_corner_indices() assumes the X/Y vectors are uniform for all rows/columns, which isn't the case when the space pixels are removed.