Differences in resampling coverage when area crosses 180 deg longitude

[pnuu]

Describe the bug
It seems that there are differences in the resampling results between satpy maint/0.8 branch (non-xarray/dask) and current master branch of satpy when the target area crosses 180 deg longitude.

To Reproduce
Resample to area definition (example below) with satpy maint/0.8 and satpy master, and compare the areas covered. The version of pyresample stays the same.

areas.def:

REGION: EPSG4326 {
   NAME:    Global equal latitude/longitude grid for global sphere
   PCS_ID:  EPSG:4326
   PCS_DEF: init=EPSG:4326
   XSIZE: 8192
   YSIZE: 4096
   AREA_EXTENT:  (-180.0, -90.0, 180.0, 90.0)
};

Test script using GOES-15 data. HIMAWARI-8/AHI data shows similar behaviour:

from satpy import Scene
import glob
fnames = glob.glob("/home/lahtinep/data/satellite/geo/L*20180806*")
glbl = Scene(reader='hrit_goes', filenames=fnames)
glbl.load([10.7])
lcl = glbl.resample('EPSG4326', resampler='nearest')
lcl.save_dataset(10.7, '/tmp/test.png')

Expected behavior
IR 10.7 um processed with maint/0.8 branch:

Actual results
Same data processed with current master branch:

Environment Info:

• SatPy Version: current master branch (0.9.1a0.dev0)
• PyResample Version: 1.10.1

[pnuu]

If the area reduction in https://github.com/pytroll/satpy/blob/master/satpy/scene.py#L897-L906 is commented out, the results are as expected.

[pnuu]

Seems to be a bug in pyresample.spherical_geometry.intersection_polygon():

In green: EPSG:4326 boundary, in blue: GOES-15 boundary, in red: intersection. The minimum longitude of the intersection is -187.376... degrees, and should be on the other edge.

[djhoese]

@pnuu Has there been any work on this on the pyresample side of things? Are resampling results the same between 0.8 and 0.9 for non-180 data?

[pnuu]

I know @mraspaud started something, but I guess he got distracted on something else. I'll check with couple different scenes tomorrow morning to check if there are differences in resampling results when the area doesn't wrap around.

[mraspaud]

I got distracted indeed, but I haven't started anything either. The point was that the current boolean union/intersection algorithm doesn't expect multiple polygons as a result, so it just returns the first intersection/union polygon it finds. The resolution of this is thus in pyresample, with possibly need for changing the api to allow returning multiple polygons to satpy, which in turn would have to handle this properly.

[ghiggi]

I finally just found also the root cause of this.
In pyresample.spherical.arc.intersections there is the need to copy self and other_arc to avoid possible external changes of such Arc(s) when the if <condition> are True ... which lead to the problem above described.
PR on spherical.py coming soon :=)

[djhoese]

🤦 the Arc properties are modified inplace

[ghiggi]

Yep :-P ... and took me 3 hours to find it out ahaha

[djhoese]

I don't think you need to copy the entire Arc, just keep a separate variable/handle for each value used later in the method.

[INElutTabile]

Hi everyone,
I'm Fabrizio Borgia from EUMETSAT, where I work at the EUMETView service. We're in the process of switching over to Satpy for generating most of our imagery, including the Geostationary ring layers (like this one).

We've run into the issue described here while working with data that crosses the 180° meridian. I noticed this issue has been merged, but we're still experiencing the problem in the version of SatPy we're using (0.43.0, should be the latest).

Any updates on this would be greatly appreciated.

Cheers and thanks!

[djhoese]

@INElutTabile Could you provide a minimal reproducible example (assuming we have equivalent data to what you're using)? Or is your example equivalent to @pnuu's original code in the first comment of this issue?

[INElutTabile]

@djhoese here is a reproducible example, and the images we produced with that, showhing the issue.

#!/usr/bin/env python

from satpy import Scene
import glob

fnames = glob.glob("data/IMG_DK01???_20231010110?_???.bz2")
scn = Scene(reader='ahi_hrit', filenames=fnames)
scn.load([10.4])
resampled = scn.resample('worldeqc3km73', resampler='nearest')
resampled.save_dataset(10.4, 'himawari_ir104_202310101100.png')

fnames = glob.glob("data/*G18*s20232831100*")
scn = Scene(reader='abi_l1b', filenames=fnames)
scn.load([10.35])
resampled = scn.resample('worldeqc3km73', resampler='nearest')
resampled.save_dataset(10.35, 'goes_ir1035_202310101100.png')

GOES

HIMAWARI

[pnuu]

@INElutTabile a workaround is to add recude_data=False to the scn.resampl() calls.

[djhoese]

Just to be clear, the himawari one is too short (the image is larger than this) and the GOES one is fine, but perhaps the default radius of influence (dynamically determined) is too low so it has large artifacts on the left. Right?

[INElutTabile]

Thanks @pnuu and @djhoese for your help - the reduce_data=False workaround in the scn.resample() call did the trick!
About the Moiré pattern like effects on the sides, that is indeed related to the radius of influence. That part is however not visible in the final composite, because it is covered or cut off in the Geostationary Ring due to the overlapping of the satellite full disks.

The updated script looks like this.

#!/usr/bin/env python

from satpy import Scene
import glob

fnames = glob.glob("data/IMG_DK01???_20231010110?_???.bz2")
scn = Scene(reader='ahi_hrit', filenames=fnames)
scn.load([0.64])
resampled = scn.resample('worldeqc3km73', resampler='nearest', reduce_data=False)
resampled.save_dataset(0.64, 'himawari_visi006_202310101100.png')

fnames = glob.glob("data/*G18*s20232831100*")
scn = Scene(reader='abi_l1b', filenames=fnames)
scn.load([13.30])
resampled = scn.resample('worldeqc3km73', resampler='nearest', reduce_data=False)
resampled.save_dataset(13.30, 'goes_ir133_202310101100.png')

Thanks again!

[arcanerr]

From #383 (comment):

Just to be clear, the himawari one is too short (the image is larger than this) and the GOES one is fine, but perhaps the default radius of influence (dynamically determined) is too low so it has large artifacts on the left. Right?

Himawari and GOES were mixed up and after you wrote your comment, @djhoese, the captions in @INElutTabile's previous comment (#383 (comment)) were swapped to fix this. So just to reduce confusion 😉: the left part of GOES is too short and Himawari has the "bite" artifact in the left part.

Both seem to be effects of the problem at hand: with reduce_data=False (thanks for the hint, @pnuu) they vanish. When reading the documentation for reduce_data (https://satpy.readthedocs.io/en/latest/api/satpy.scene.html#satpy.scene.Scene.resample) it starts to make sense, why wrongly calculating an intersection polygon causes this issue.

Personally, I don't like that reduce_data=True is the default setting, but in any case it seems advisable to set reduce_data=False if the target area covers (almost) the entire earth or the source area is known to be completely contained in the (inverse projected) target area -- and/or if the source data overlaps 180 degree longitude (that is the bug).

Either way, the problem persists and could also occur with other target projections/data that cover the 180 degree longitude. Might be worth to test, if it also happens at other longitudes than 180 deg, if the target projection is cut at a different longitude.

[djhoese]

Thanks for the clarification and summary @arcanerr. The problem with reduce_data is overall a bug in the intersection/boundary polygon code as you've pointed out. In my opinion reduce_data=False as the default is not a good option. It would result in much worse performance for anyone doing any type of resampling to a smaller than full disk area (I haven't proven this in a long time, but theoretically it should be much worse). Additionally, I make the assumption that most people who are doing "full image" images are using native resampling which should not suffer/use reduce_data at all. It seems to be the biggest problem for large target areas and/or target areas that are near the edge of the disk (which for people like @pnuu is most areas).