pytroll / pyresample

Geospatial image resampling in Python
http://pyresample.readthedocs.org
GNU Lesser General Public License v3.0
346 stars 94 forks source link

Resampling swath to another swath: no attribute get_proj_coords #436

Open adriaat opened 2 years ago

adriaat commented 2 years ago

According to the Resampling of swath data docs:

Pyresample can be used to resample a swath dataset to a grid, a grid to a swath or a swath to another swath.

but if

from pyresample import geometry
from pyresample.bilinear import XArrayBilinearResampler

# Missing variable definitions omitted for simplicity,
# but all 2D arrays, and radius_of_influence a number > 0 representing metres
source_def = geometry.SwathDefinition(lons=lons_s, lats=lats_s)
target_def = geometry.SwathDefinition(lons=lons_t, lats=lats_t)
resampler = XArrayBilinearResampler(source_def, target_def, radius_of_influence)
result = resampler.resample(data)

then the following error is raised:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/home/amell/test/whiteboard.ipynb Cell 27' in <cell line: 1>()
----> [1](vscode-notebook-cell://home/amell/test/whiteboard.ipynb#ch0000027vscode-remote?line=0) resampled_surface_precipitation = resampler.resample(ds_gpm_subset.surface_precipitation.values)

File ~/miniconda3/envs/test/lib/python3.10/site-packages/pyresample/bilinear/xarr.py:71, in XArrayBilinearResampler.resample(***failed resolving arguments***)
     69 """Resample the given data."""
     70 del nprocs
---> 71 self.get_bil_info()
     72 return self.get_sample_from_bil_info(data, fill_value=fill_value, output_shape=None)

File ~/miniconda3/envs/test/lib/python3.10/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/envs/test/lib/python3.10/site-packages/pyresample/bilinear/_base.py:165, in BilinearBase._get_fractional_distances(self)
    164 def _get_fractional_distances(self):
--> 165     out_x, out_y = self._get_output_xy()
    166     # Get the four closest corner points around each output location
    167     corner_points, self._index_array = \
    168         _get_four_closest_corners(*self._get_input_xy(),
    169                                   out_x, out_y,
...
    217 def _get_output_xy(target_geo_def):
--> 218     out_x, out_y = target_geo_def.get_proj_coords(chunks=CHUNK_SIZE)
    219     return da.compute(np.ravel(out_x), np.ravel(out_y))

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

The same happens if NumpyBilinearResampler is used in place of XArrayBilinearResampler.

I expected pyresample to be consistent with the docs. On the other hand, in this example I expect that radius_of_influence is used on the WSG84, and not on some Cartesian coordinates.

Versions used:

djhoese commented 2 years ago

That sentence describes pyresample as a whole which does include resampling methods for swath to swath resampling, but not every resampling algorithm implemented in pyresample. The relatively new bilinear resampling classes do not currently support using a swath definition as the target area. I will agree that the documentation is terrible which is why I'm rewriting it in #434.

On the other hand, in this example I expect that radius_of_influence is used on the WSG84, and not on some Cartesian coordinates.

I understand how this could be confusing and we've tried updating documentation to reflect what radius of influence actually is. The units/coordinates are geocentric cartesian (euclidean) coordinate because that's what the internal KDTree structure uses. This allows us to avoid issues with pixels being over the edge of the input or output projections or hitting the anti-meridian of a projection.

jdldeauna commented 2 years ago

Hi @djhoese,

I was just wondering if the bilinear resampling classes also do not support a GridDefinition as target data. When I use code similar to one posted above:

source_def = geometry.SwathDefinition(lons=lons_s, lats=lats_s)
target_def = geometry.GridDefinition(lons=new_grid_lons, lats=new_grid_lats)
resampler = XArrayBilinearResampler(source_def, target_def, radius_of_influence)
result = resampler.resample(data)

results in the ff error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [17], in <cell line: 1>()
----> 1 result = resampler.resample(da)

File /srv/conda/envs/notebook/lib/python3.9/site-packages/pyresample/bilinear/xarr.py:71, in XArrayBilinearResampler.resample(***failed resolving arguments***)
     69 """Resample the given data."""
     70 del nprocs
---> 71 self.get_bil_info()
     72 return self.get_sample_from_bil_info(data, fill_value=fill_value, output_shape=None)

File /srv/conda/envs/notebook/lib/python3.9/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 /srv/conda/envs/notebook/lib/python3.9/site-packages/pyresample/bilinear/_base.py:165, in BilinearBase._get_fractional_distances(self)
    164 def _get_fractional_distances(self):
--> 165     out_x, out_y = self._get_output_xy()
    166     # Get the four closest corner points around each output location
    167     corner_points, self._index_array = \
    168         _get_four_closest_corners(*self._get_input_xy(),
    169                                   out_x, out_y,
    170                                   self._neighbours, self._index_array)

File /srv/conda/envs/notebook/lib/python3.9/site-packages/pyresample/bilinear/xarr.py:91, in XArrayBilinearResampler._get_output_xy(self)
     90 def _get_output_xy(self):
---> 91     out_x, out_y = _get_output_xy(self._target_geo_def)
     92     out_x = out_x[self._valid_output_indices]
     93     out_y = out_y[self._valid_output_indices]

File /srv/conda/envs/notebook/lib/python3.9/site-packages/pyresample/bilinear/xarr.py:218, in _get_output_xy(target_geo_def)
    217 def _get_output_xy(target_geo_def):
--> 218     out_x, out_y = target_geo_def.get_proj_coords(chunks=CHUNK_SIZE)
    219     return da.compute(np.ravel(out_x), np.ravel(out_y))

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

Versions used:

djhoese commented 2 years ago

@jdldeauna That is correct that it doesn't. However, making it support a GridDefinition is much easier than a SwathDefinition. Swaths by definitions are not uniformly spaced. A Grid is by definition uniformly spaced. So you could:

  1. Submit a pull request to pyresample to add the necessary method(s) to GridDefinition to return slices of the original lon/lat grids as needed by the resampler.
  2. Change to creating an AreaDefinition instead of a GridDefinition. If you define your AreaDefinition with a lon/lat CRS like "EPSG:4326" as your "proj_dict". Make sure your extents are calculated to be the outer edges of your lon/lat array coordinates (which are defined as the centers of the pixels). In my opinion and as I think of where I want pyresample 2.0 to go, I think GridDefinition is just a convenience method/class around the AreaDefinition class where it should just be doing these calculations for you.