Geographic EWA projection for swaths crossing the anti-meridian omits values for 90 ≤ longitude (degrees east) ≤ 180.

[owenlittlejohns]

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

Input data swath crossing anti-meridian: Available here, also attached to this issue.

Note, this summarises code from various parts of a larger service - so is a little rough, particularly with handling the fill value.

from pyresample import AreaDefinition, SwathDefinition
from pyresample.ewa import fornav, ll2cr
from pyresample.utils import check_and_wrap
from netCDF4 import Dataset
import numpy as np

INPUT_FILE_PATH = '20210817140001-JPL-L2P_GHRSST-SSTskin-MODIS_A-N-v02.0-fv01.0.nc'
TARGET_DIMENSIONS = (433, 7265)
TARGET_EXTENTS = (-179.992447, 31.104006, 179.997894, 52.569942)

with Dataset(INPUT_FILE_PATH, 'r') as input_dataset:
    # Get fill value
    unscaled_fill_value = variable.getncattr('_FillValue')
    scale_factor = variable.getncattr('scale_factor')
    add_offset = variable.getncattr('add_offset')
    variable_fill_value = (unscaled_fill_value * scale_factor) + add_offset

    variable_data = input_dataset['sea_surface_temperature'][0][:].filled(
        fill_value=variable_fill_value
    )
    input_swath = SwathDefinition(lats=wrapped_lats, lons=wrapped_lons)
    target_grid = AreaDefinition.from_extent('lat, lon', 'EPSG:4326',
                                             TARGET_DIMENSIONS,
                                             TARGET_EXTENTS)

    ewa_info = ll2cr(input_swath, target_grid)
    _, results = fornav(ewa_info[1], ewa_info[2], target_grid,
                        variable_data, maximum_weight_mode=False)

    np.nan_to_num(results, nan=variable_fill_value, copy=False)

Problem description

When taking an input swath that crosses the anti-meridian and reformatting it to a whole-Earth geographic grid, the EWA algorithm does not include data where 90 < longitude (degrees east) < 180. (We've confirmed that this behaviour is limited to the EWA algorithm, other interpolation methods work nicely).

My colleagues and I have looked into this a bit, and think the issue may be in ll2cr_static. There's a check for anti-meridian crossing. If the pixel is determined to cross the anti-meridian, the x coordinate (longitude in this case) is incremented by the projection circumference (360 degrees for geographic):

elif proj_circum != 0 and abs(x_tmp - origin_x) >= (0.75 * proj_circum):
    # if x is more than 75% around the projection space, it is probably crossing the anti-meridian
    x_tmp += proj_circum

In our case, we have a whole-Earth output grid, so origin_x = -180 and proj_circum = 360.0. The anti-meridian crossing condition for these values is satisfied for longitude ≥ 90 (degrees east). In the case of longitude = 90 degrees east, this shift means longitude_incremented = 450 degrees east. This is likely leading to the affected longitudes being considered outside of the target grid.

Our team was wondering if this could be mitigated by an additional check to the elif condition that ensures the x coordinate is less than the lefthand edge of the target grid: x_tmp < origin_x. Alternatively, maybe the projection could be checked to see if the target grid is geographic? (Or maybe a combination of both checks)

Expected Output

This is the input used in the sample code above. The output should resemble the input, but be geographically gridded to the target area:

Actual Result

The output below shows the non-fill pixels once the sea_surface_temperature has been projected. The output target area is essentially whole-Earth, but the pixels west of the anti-meridian are filled.

Panoply georeferenced image	Gridded contour plot using target grid

Versions of Python, package at hand and relevant dependencies

Tested under pyresample v1.17.0 and v1.20.0 in Python 3.7.

[djhoese]

Thanks for filing the issue!

It has probably been nearly 10 years since I originally wrote this code in my Polar2Grid project so you'll have to forgive my lack of memory for why it is the way it is. My gut reaction is that this if statement and the logic for projection circumference just get removed. Do you think you and your colleagues could try installing a development version of pyresample with that elif commented out and see if it fixes the issue? That would be the first step in moving forward with this.

The two things I can think of for why this is here is:

1. Avoid data being interpolated from the left edge of the output grid to the right edge of the output grid. I don't think this is possible with the EWA algorithm, but it is hard to remember.
2. Avoid data that is just over the antimeridian not being used for interpolation for the output pixels just on the other edge of the antimeridian. They are geographically next to each other so it would be nice if they influenced each other. I don't think this can be dealt with though.

I just need to make sure I walk through the various use cases of doing things like a global mercator area versus global lon/lat area versus a lon/lat area that has -180/180 as its origin/center.

[owenlittlejohns]

Thanks for the quick response. We'll try installing a development version commenting out the elif and see what we get! (Not sure on the timeframe, but will keep you posted)

[owenlittlejohns]

@djhoese - I just got chance to try out a development version of pyresample locally. I tried two things:

• Removing the elif, as suggested.
• Combining the second clause of the elif statement with a check for x_tmp < origin_x. This made the full statement:

elif proj_circum != 0 and (origin_x - x_tmp) >= (0.75 * proj_circum):
    x_tmp += proj_circum

I tried both of these approaches with my example file (as shown above) for some of the target grids we see fairly often:

• Geographic ("EPGS:4326").
• Cylindrical Equal Area ("EPSG:6933").
• Lambert Azimuthal Equal Area ("+proj=laea").
• Mercator ("EPSG:3857").

In all cases, and both approaches, the output contained the data on both sides of the anti-meridian as expected!

Panoply Georeferenced Image	Gridded contour plot using target grid

[djhoese]

Great! How would you feel about making a pull request with the first change (removing the whole elif)?

[owenlittlejohns]

Sounds good to me!