Better handling of anti-meridian cases

[djhoese]

Sorry for not following the template, but I'm not sure this fits.

Currently the DynamicAreaDefinition.freeze method has an antimeridian_mode keyword argument that I added for handling cases where the provided lon/lats cross the anti-meridian (-180/180 in lonlat space). This kwarg can be:

1 .modify_crs (default): Add a +pm=180 to shift the prime meridian of the projection to 180 degrees. This means the extents surround 0 instead of being discontinuous over -180/180. See downsides below.
2. modify_extents: Mess with extents of lon/lat projections so the longitudes are 0-360 instead of a -180-180 range. See downsides below.
3. global_extents: Change the extents to be the entire globe. Theoretically if resampling handling this properly you get a full globe with the input data split into two pieces.

The modify_crs is problematic because many tools don't respect +pm=180 (see https://proj.org/en/9.5/usage/projections.html#prime-meridian). The modify_extents is a problem because it seems some (all?) of our resamplers can't handle longitudes outside the -180-180 range either deliberately filtering them as invalid (kdtree nearest neighbor does this) or suffering from PROJ transformation automatically wrapping longitudes to the -180-180 range.

For these reasons I propose a few important changes:

1. Deprecate modify_exents and combine it with modify_crs, but at the same time change modify_crs to use +lon_wrap=180 instead of +pm=180. See https://proj.org/en/9.5/usage/projections.html#longitude-wrapping which I think is what we normally want anyway.
2. Change kdtree resamplers to only filter lon/lats >=1e30 instead of doing the full -180/180 and -90/90 validation. This would allow PROJ to handle the wrapping.
3. Update EWA resampler to use pyproj's Transformer class and be more explicit about convert lon/lats to the target projection. This should hopefully gracefully use the +lon_wrap=180 which I don't think is working right now.

Questions:

• What does get_lonlats for an area return that has +lon_wrap=180. Should it return 0-360 or -180/180 even though they wouldn't be contiguous? In the past I think we updated this method to assert that lon/lats are always -180/180.
• Edit: Does +pm=180 get us anything (is it useful to us or our users) that +lon_wrap=180 doesn't?

[sfinkens]

I'm not sure I understand all aspects of the problem. But from my limited view I would expect get_lonlats to return 0-360 when I set +wrap_lon=180 in the projection definition.

I made a quick test with an area def that goes from 120 to -120 degrees longitude.

import pyresample
import numpy as np

def test_wrap(wrap_key, extent):
    area_global = pyresample.create_area_def(
        area_id="test",
        projection={'datum': 'WGS84', 'proj': 'longlat'},
        area_extent=[-180, -60, 180, 60],
        resolution=60
    )
    area_wrapped = pyresample.create_area_def(
        area_id="test",
        projection={'datum': 'WGS84', wrap_key: 180.0, 'proj': 'longlat'},
        area_extent=extent,
        resolution=60
    )

    # Compute coordinates
    lons, lats = area_wrapped.get_lonlats()
    x, y = area_wrapped.get_proj_vectors()

    # Resample fake data
    data = np.ones(area_wrapped.shape)
    resampled = pyresample.kd_tree.resample_nearest(area_wrapped, data, area_global, radius_of_influence=50000)

    # Print results
    print(lons)
    print(x)
    print(resampled)

With pm:

>>> test_wrap("pm", [-60, -60, 60, 60])
[[ 150. -150.]
 [ 150. -150.]]
[-30.  30.]
[[1. 0. 0. 0. 0. 1.]
 [1. 0. 0. 0. 0. 1.]]

With lon_wrap (not sure if I got the extent right):

test_wrap("lon_wrap", [120, 30, 240, 150])
[[150. 210.]
 [150. 210.]]
[150. 210.]
[[0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]]

So the resampling obviously doesn't work. But what I like about lon_wrap is that lon and x are identical. That difference often confused me with pm.

[djhoese]

Yeah, pm can be considered a "workaround" at best. You're literally changing the projection's "home" and users then can't expect this to be the same as what they put in (projection in != projection out).

[ghiggi]

Hi @djhoese ! I just want to point out a few details I have in my notes that might help clarifying/solving this issue.

• Currently, the k-d tree is built using geocentric x, y, z coordinates. In the kdtree module, we first extract lon,lat coordinates with get_lonlat for each source and destination area, then convert them to x, y, z.

• The lon,lat coordinates are currently required to retrieve x, y, z coordinates and for performing data reduction when reduce_data=True.

• AreaDefinition.get_lonlat returns lon,lat in the geodetic CRS with +pm=0

• SwathDefinition.get_lonlat returns the original lon,lat of the source CRS, which might not be +pm=0. This difference cause bad results for a SwathDefinition if the CRS prime-meridian is not zero. This because the current conversion to x, y, z in the kdtree module assumes lon,lat with a prime meridian of 0 and values ranging from -180 to 180 degrees (see https://github.com/pytroll/pyresample/blob/main/pyresample/future/resamplers/nearest.py#L61, https://github.com/pytroll/pyresample/blob/main/pyresample/future/resamplers/_transform_utils.py#L22, https://github.com/pytroll/pyresample/blob/main/pyresample/kd_tree.py#L533)

• To avoid multiple coordinate conversions in the kdtree module and remove duplicated and fragile code, we could add an xyz = area.get_geocentric_coordinates(data_slice=None, chunks=None) method to the area geometries. This method would use the dask map_blocked pyproj.Transform and could be called directly in kdtree.

• If get_geocentric_coordinates would be used, computations of lon,lat would be required only if reduce_data=True.

• However, reduce_data=True logic requires lon/lat values between -180 and 180 degrees and a contiguous longitude array. Without modifying the logic in data_reduce.py reduce_data=True will continue to yield incorrect results if an area crosses the antimeridian.

• Data reduction might be made more efficient by using the source_area.get_area_slices(destination_area) as currently done in Satpy (see this code section).

• Satpy seems to not use anymore the reduce_data code of pyresample, although some apparent old unused reduce_data argument are still around (see here).

• Using source_area.get_area_slices(destination_area) would utilize the boundary logic and spherical geometry operations, avoiding the need for full lon,lat arrays computations.

[TomLav]

Hello. I admit I do not know exactly what this issue / enhancement is about. But since the antimeridian and data_reduce are mentioned, I would link to another github issue about these two topics from back in 2017.
#61

Maybe the topics are un-related, maybe the issue from 2017 should have been closed. But I know that I from time to time get unexpected behaviour with kdtree and grids having the pole / antimeridian. And that setting data_reduce to False is my quick fix.

Sorry if I am off-topic here.