AttributeError when trying to resample grid data to swath data using NumpyBilinearResampler.

[danschef]

Code Sample, a minimal, complete, and verifiable piece of code

import numpy as np
from pyresample.bilinear import NumpyBilinearResampler
from pyresample import geometry

area_def = \
    geometry.AreaDefinition('areaD',
                            'Europe (3km, HRV, VTC)',
                            'areaD',
                            {'a': '6378144.0', 'b': '6356759.0',
                             'lat_0': '50.00', 'lat_ts': '50.00',
                             'lon_0': '8.00', 'proj': 'stere'},
                            800, 800,
                            [-1370912.72, -909968.64,
                             1029087.28, 1490031.36])
swath_data = np.fromfunction(lambda y, x: y * x, (500, 100))
lons = np.fromfunction(lambda y, x: 3 + x * 0.1, (500, 100))
lats = np.fromfunction(lambda y, x: 75 - y * 0.1, (500, 100))
swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
grid_data_resamp = NumpyBilinearResampler(swath_def, area_def, 30e3).resample(swath_data)
swath_data_resamp = NumpyBilinearResampler(area_def, swath_def, 30e3).resample(grid_data_resamp)  ## <- fails

Problem description

I am using pyresample for transforming swath data to grid data and vice versa. Using gauss resampling or nearest neighbor resampling, everything works fine in both directions. The code example for bilinear resampling in the docs is also working without any issues. However, it only contains the case of resampling a swath array into a grid array and not vice versa. But I also need the other direction.

When I swap the source and target geometry definition which are passed to NumpyBilinearResampler (as in the code sample above), pyresample raises the AttributError below. I used the example code snippet from the docs, I just added the resampling in the other direction (last line, which fails).

Am I doing something wrong here? If yes, could someone show me how to use NumpyBilinearResampler to resample grid data to swath data?

Expected Output

I would expect that the bilinear resampling produces a similar result like the gauss resampling where it is easiliy possible to resample from swath to grid data and vice versa (code is taken from the docs, only variable names are changed):

import numpy as np
from pyresample import kd_tree, geometry
area_def = \
    geometry.AreaDefinition(
        'areaD',
        'Europe (3km, HRV, VTC)',
        'areaD',
        {'a': '6378144.0', 'b': '6356759.0',
         'lat_0': '50.00', 'lat_ts': '50.00',
         'lon_0': '8.00', 'proj': 'stere'},
        800, 800,
        [-1370912.72, -909968.64,
         1029087.28, 1490031.36])
swath_data = np.fromfunction(lambda y, x: y * x, (50, 10))
lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
grid_data_resamp = kd_tree.resample_gauss(swath_def, swath_data, area_def,
                                          radius_of_influence=50000, sigmas=25000)
swath_data_resamp = kd_tree.resample_gauss(area_def, grid_data_resamp, swath_def,
                                           radius_of_influence=50000, sigmas=25000)

When plotted, the arrays look like this: Note, that swath_data_resamp has nearly the same values like the original swath_data array within the extent given to geometry.AreaDefinition. This is exactly how the result should look like.

Actual Result, Traceback if applicable

NumpyBilinearResampler.resample() raises AttributeError: 'SwathDefinition' object has no attribute 'get_proj_coords':

/home/gfz-fe/miniconda3/envs/pyresample117/lib/python3.9/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
  proj_string = self.to_proj4()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-41-e1c58470d841> in <module>
     18 swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     19 grid_data_resamp = NumpyBilinearResampler(swath_def, area_def, 30e3).resample(swath_data)
---> 20 swath_data_resamp = NumpyBilinearResampler(area_def, swath_def, 30e3).resample(grid_data_resamp)

~/miniconda3/envs/pyresample117/lib/python3.9/site-packages/pyresample/bilinear/_numpy_resampler.py in resample(self, data, fill_value, nprocs)
    247         else:
    248             kdtree_class = KDTree
--> 249         self.get_bil_info(kdtree_class=kdtree_class, nprocs=nprocs)
    250         return self.get_sample_from_bil_info(data, fill_value=fill_value, output_shape=None)
    251 

~/miniconda3/envs/pyresample117/lib/python3.9/site-packages/pyresample/bilinear/_base.py in get_bil_info(self, kdtree_class, nprocs)
    114 
    115         # Calculate vertical and horizontal fractional distances t and s
--> 116         self._get_fractional_distances()
    117         self._get_target_proj_vectors()
    118         self._get_slices()

~/miniconda3/envs/pyresample117/lib/python3.9/site-packages/pyresample/bilinear/_base.py in _get_fractional_distances(self)
    159 
    160     def _get_fractional_distances(self):
--> 161         out_x, out_y = self._get_output_xy()
    162         # Get the four closest corner points around each output location
    163         corner_points, self._index_array = \

~/miniconda3/envs/pyresample117/lib/python3.9/site-packages/pyresample/bilinear/_base.py in _get_output_xy(self)
    169 
    170     def _get_output_xy(self):
--> 171         return _get_output_xy(self._target_geo_def)
    172 
    173     def _get_input_xy(self):

~/miniconda3/envs/pyresample117/lib/python3.9/site-packages/pyresample/bilinear/_base.py in _get_output_xy(target_geo_def)
    264 
    265 def _get_output_xy(target_geo_def):
--> 266     out_x, out_y = target_geo_def.get_proj_coords()
    267     return np.ravel(out_x),  np.ravel(out_y)
    268 

AttributeError: 'SwathDefinition' object has no attribute 'get_proj_coords'

Arrays look like this (swath_data_resamp is missing):

Versions of Python, package at hand and relevant dependencies

pyyresample 1.17.0, installed via conda install -c conda-forge pyresample; Python 3.9

[pnuu]

Hi, and sorry for the late reply.

This is an unfortunate side-effect of how the algorithm we use for bilinear interpolation works. For the bilinear interpolation we need to find four closest data points from the source data that surround the requested target location. The distance the algorithm uses is in the target area's projection coordinates. Swath data doesn't have a valid projection, thus the target area's projection coordinates can't be calculated.

I think it could be possible to adapt the resampler to use 3D geocentric cartesian (X, Y, Z) coordinates. But at the least this shortcoming should be indicated in the documentation until we can address the issue.

[danschef]

Thanks for the reply. I agree, if there is a limitation of the current bilinear algorithm, this should be mentioned in the docs. Additionally, I would also return a meaningful exception in case NumpyBilinearResampler is called to resample grid data to swath data.

Honestly, I don´t really get why you need a valid projection for the swath data. I mean with swath_def = geometry.SwathDefinition(lons=lons, lats=lats) you have a clear mapping between cartesian image coordinates and longitude/latitude coordinates. If needed, you can also easily transform these lon/lat coordinates to the projection of the grid data to have both in the same projection. What makes it difficult to determine the 4 closest pixels within the grid data array for a given lon/lat coordinate?

[danschef]

@pnuu So far, I worked around this issue by pinning pyresample to <1.17. This now increasingly causes trouble and blocks me from building a subsequent package for Python 3.10.

What do you think is the effort to adapt the bilinear resampler to use cartesian coordinates instead of target area's projection coordinates when resampling grid data to swath data? Or is there any workaround I could use for now to make this work?