Repository: pytroll/pyresample (language: python)

Issue:
Geographic EWA projection for swaths crossing the anti-meridian omits values for 90 ≤ longitude (degrees east) ≤ 180.
#### Code Sample, a minimal, complete, and verifiable piece of code

Input data swath crossing anti-meridian: [Available here](https://podaac-tools.jpl.nasa.gov/drive/files/allData/ghrsst/data/GDS2/L2P/MODIS_A/JPL/v2019.0/2021/229/20210817140001-JPL-L2P_GHRSST-SSTskin-MODIS_A-N-v02.0-fv01.0.nc), also [attached to this issue](https://github.com/pytroll/pyresample/files/7086053/20210817140001-JPL-L2P_GHRSST-SSTskin-MODIS_A-N-v02.0-fv01.0.nc.gz).

Note, this summarises code from various parts of a larger service - so is a little rough, particularly with handling the fill value.

```python
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`](https://github.com/pytroll/pyresample/blob/main/pyresample/ewa/_ll2cr.pyx#L255). 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):

```python
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:

<img src="https://user-images.githubusercontent.com/7788154/131731132-49553619-1170-4eae-bbf1-8e67b5d92db4.png" width="400px">


#### 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 |
|---|---|
| <img src="https://user-images.githubusercontent.com/7788154/131731274-562ed401-6439-40f5-a2fa-444be758ac46.png" width="350px"> | <img src="https://user-images.githubusercontent.com/7788154/131731288-13b63970-4426-438c-bd47-cd6829b38927.png" width="350px"> |


#### Versions of Python, package at hand and relevant dependencies

Tested under `pyresample` v1.17.0 and v1.20.0 in Python 3.7.


Discussion / hints:
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.
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)
@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:

```python
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 |
| --- | --- |
| <img src="https://user-images.githubusercontent.com/7788154/132911011-29943801-2b32-4d58-a127-8c4e5bef127d.png" width="350px"> | <img src="https://user-images.githubusercontent.com/7788154/132911014-49cdce92-289e-41ca-80a4-8fb802b063db.png" width="350px"> |
Great! How would you feel about making a pull request with the first change (removing the whole `elif`)?
Sounds good to me!

Fix patch:
diff --git a/pyresample/ewa/_ll2cr.pyx b/pyresample/ewa/_ll2cr.pyx
index d009010e6..4a42cab1f 100644
--- a/pyresample/ewa/_ll2cr.pyx
+++ b/pyresample/ewa/_ll2cr.pyx
@@ -232,7 +232,6 @@ def ll2cr_static(numpy.ndarray[cr_dtype, ndim=2] lon_arr, numpy.ndarray[cr_dtype
     cdef tuple projected_tuple = p(lon_arr, lat_arr)
     cdef cr_dtype [:, ::1] rows_out = projected_tuple[1]
     cdef cr_dtype [:, ::1] cols_out = projected_tuple[0]
-    cdef double proj_circum = projection_circumference(p)
 
     # indexes
     cdef unsigned int row
@@ -252,9 +251,6 @@ def ll2cr_static(numpy.ndarray[cr_dtype, ndim=2] lon_arr, numpy.ndarray[cr_dtype
                 lon_arr[row, col] = fill_in
                 lat_arr[row, col] = fill_in
                 continue
-            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
 
             x_tmp = (x_tmp - origin_x) / cell_width
             y_tmp = (y_tmp - origin_y) / cell_height


Test patch:
diff --git a/pyresample/test/test_ewa_ll2cr.py b/pyresample/test/test_ewa_ll2cr.py
index f3e666ae6..2fb7f836e 100644
--- a/pyresample/test/test_ewa_ll2cr.py
+++ b/pyresample/test/test_ewa_ll2cr.py
@@ -50,6 +50,17 @@
     "proj4_definition": "+proj=lcc +a=6371200 +b=6371200 +lat_0=25 +lat_1=25 +lon_0=-95 +units=m +no_defs",
 }
 
+static_geo_whole_earth = {
+    "grid_name": "test_geo",
+    "origin_x": -180.0,
+    "origin_y": 40,
+    "width": 361,
+    "height": 181,
+    "cell_width": 1.0,
+    "cell_height": 1.0,
+    "proj4_definition": "+proj=longlat +datum=WGS84 +no_defs +type=crs",
+}
+
 
 class TestLL2CRStatic(unittest.TestCase):
     def test_lcc_basic1(self):
@@ -66,9 +77,32 @@ def test_lcc_basic1(self):
         w = grid_info["width"]
         h = grid_info["height"]
         points_in_grid = _ll2cr.ll2cr_static(lon_arr, lat_arr, fill_in, proj_str,
-                                                               cw, ch, w, h, ox, oy)
+                                             cw, ch, w, h, ox, oy)
         self.assertEqual(points_in_grid, lon_arr.size, "all these test points should fall in this grid")
 
+    def test_geo_antimeridian(self):
+        """ Ensure output for anti-meridian crossing input includes all points. """
+        from pyresample.ewa import _ll2cr
+        lon_arr = create_test_longitude(160.0, 200.0, (50, 100), dtype=np.float64)
+        lat_arr = create_test_latitude(40.0, 60.0, (50, 100), dtype=np.float64)
+
+        # wrap values in lon_arr so -180 ≤ longitude (degrees east) ≤ 180
+        lon_arr[lon_arr > 180] -= 360.0
+
+        grid_info = static_geo_whole_earth.copy()
+        fill_in = np.nan
+        proj_str = grid_info['proj4_definition']
+        cw = grid_info['cell_width']
+        ch = grid_info['cell_height']
+        ox = grid_info['origin_x']
+        oy = grid_info['origin_y']
+        w = grid_info['width']
+        h = grid_info['height']
+        points_in_grid = _ll2cr.ll2cr_static(lon_arr, lat_arr, fill_in, proj_str,
+                                             cw, ch, w, h, ox, oy)
+        self.assertEqual(points_in_grid, lon_arr.size,
+                         'all these test points should fall in this grid')
+
     def test_lcc_fail1(self):
         from pyresample.ewa import _ll2cr
         lon_arr = create_test_longitude(-15.0, 15.0, (50, 100), dtype=np.float64)