Resampling with gradient search results in interlacing image artefacts if corner pixels have no data

[gerritholl]

Describe the bug

When resampling with gradient search to my target area, the resulting image exhibits interlacing image artefacts. In gwenview, the image is shown with large slices, although the exact slicing depends on the zoom level. In the gimp, the image is shown interlaced with black lines (mode L) or semitransparent (mode LA).

When I use nearest neighbour resampling, the resulting image looks correct.

To Reproduce

import hdf5plugin
from satpy import Scene
from satpy.utils import debug_on; debug_on()
from glob import glob
from pyresample import create_area_def
fci_files = glob("/media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928*_C_0076_*.nc")
ar = create_area_def(
        "test", 4087,
        area_extent=[-10000000, -10000000, 10000000, 10000000],
        resolution=3000)
sc = Scene(filenames={"fci_l1c_nc": fci_files})
sc.load(["vis_06"])
ls = sc.resample(ar, resampler="gradient_search")
ls.save_datasets(filename="/tmp/test/{platform_name}-{sensor}-{start_time:%Y%m%d%H%M}-{area.area_id}-{name}-test.tif", fill_value=0)

Expected behavior

I expect an image that displays correctly in standard tools.

Actual results

Selected console output

[WARNING: 2024-04-02 13:10:53 : root] shape found from radius and resolution does not contain only integers: (6666.666666666667, 6666.666666666667)
Rounding shape to (6667, 6667) and resolution from (3000.0, 3000.0) meters to (2999.850007499625, 2999.850007499625) meters
[DEBUG: 2024-04-02 13:10:53 : satpy.readers.yaml_reader] Reading ('/data/gholl/checkouts/satpy/satpy/etc/readers/fci_l1c_nc.yaml',)
<...snip...>
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.fci_l1c_nc] Reading vis_06_pixel_quality from /media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928123928_IDPFI_VAL_20230928123720_20230928123755_N__C_0076_0031.nc
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.fci_l1c_nc] Reading vis_06_pixel_quality from /media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928123939_IDPFI_VAL_20230928123729_20230928123803_N__C_0076_0032.nc
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.fci_l1c_nc] Reading vis_06_pixel_quality from /media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928123953_IDPFI_VAL_20230928123737_20230928123819_N__C_0076_0033.nc
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.fci_l1c_nc] Reading vis_06_pixel_quality from /media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928124003_IDPFI_VAL_20230928123755_20230928123835_N__C_0076_0034.nc
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.fci_l1c_nc] Reading vis_06_pixel_quality from /media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928124008_IDPFI_VAL_20230928123811_20230928123842_N__C_0076_0035.nc
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.fci_l1c_nc] Reading vis_06_pixel_quality from /media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928124015_IDPFI_VAL_20230928123819_20230928123856_N__C_0076_0036.nc
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.fci_l1c_nc] Reading vis_06_pixel_quality from /media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928124026_IDPFI_VAL_20230928123835_20230928123908_N__C_0076_0037.nc
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.fci_l1c_nc] Reading vis_06_pixel_quality from /media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928124030_IDPFI_VAL_20230928123849_20230928123914_N__C_0076_0038.nc
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.fci_l1c_nc] Reading vis_06_pixel_quality from /media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928124033_IDPFI_VAL_20230928123856_20230928123919_N__C_0076_0039.nc
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.fci_l1c_nc] Reading vis_06_pixel_quality from /media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/09/28/12/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY---NC4E_C_EUMT_20230928124035_IDPFI_VAL_20230928123908_20230928123924_N__C_0076_0040.nc
[DEBUG: 2024-04-02 13:11:07 : satpy.readers.yaml_reader] Requested orientation for Dataset None is 'native' (default). No flipping is applied.
[DEBUG: 2024-04-02 13:11:07 : satpy.scene] Resampling DataID(name='vis_06', wavelength=WavelengthRange(min=0.59, central=0.64, max=0.69, unit='µm'), resolution=1000, calibration=<1>, modifiers=())
[DEBUG: 2024-04-02 13:11:07 : satpy.scene] Data reduction disabled by the user
[DEBUG: 2024-04-02 13:11:07 : pyresample.gradient] /!\ Instantiating an experimental GradientSearch resampler /!\
[WARNING: 2024-04-02 13:11:07 : pyresample.resampler] The input area chunks are large. This usually means that the input area is of much higher resolution than the output area. You can reduce the chunks passed, and ponder whether you are using the right resampler for the job.
[WARNING: 2024-04-02 13:11:07 : pyresample.resampler] The input area chunks are large. This usually means that the input area is of much higher resolution than the output area. You can reduce the chunks passed, and ponder whether you are using the right resampler for the job.
[WARNING: 2024-04-02 13:11:07 : pyresample.resampler] The input area chunks are large. This usually means that the input area is of much higher resolution than the output area. You can reduce the chunks passed, and ponder whether you are using the right resampler for the job.
[DEBUG: 2024-04-02 13:11:07 : satpy.scene] Resampling DataID(name='vis_06_pixel_quality', resolution=1000, modifiers=())
[DEBUG: 2024-04-02 13:11:07 : satpy.scene] Data reduction disabled by the user
[DEBUG: 2024-04-02 13:11:07 : satpy.writers] Reading ['/data/gholl/checkouts/satpy/satpy/etc/writers/geotiff.yaml']
[DEBUG: 2024-04-02 13:11:07 : rasterio.session] Could not import boto3, continuing with reduced functionality.
[DEBUG: 2024-04-02 13:11:07 : satpy.writers] Adding enhancement configuration from file: /data/gholl/checkouts/satpy/satpy/etc/enhancements/generic.yaml
[DEBUG: 2024-04-02 13:11:07 : satpy.writers] Adding enhancement configuration from file: /home/gholl/Arbeit/checkouts-perforce/config/enhancements/generic.yaml
[DEBUG: 2024-04-02 13:11:07 : satpy.writers] Adding enhancement configuration from file: /home/gholl/checkouts/pytroll-dwd-config/enhancements/generic.yaml
[DEBUG: 2024-04-02 13:11:07 : satpy.writers] Data for DataID(name='vis_06', wavelength=WavelengthRange(min=0.59, central=0.64, max=0.69, unit='µm'), resolution=1000, calibration=<1>, modifiers=()) will be enhanced with options:
	[{'name': 'linear_stretch', 'method': <function stretch at 0x7f9b09cdf380>, 'kwargs': {'stretch': 'crude', 'min_stretch': 0.0, 'max_stretch': 100.0}}, {'name': 'gamma', 'method': <function gamma at 0x7f9b09cdfd80>, 'kwargs': {'gamma': 1.5}}]
[DEBUG: 2024-04-02 13:11:07 : trollimage.xrimage] Applying stretch crude with parameters {'min_stretch': 0.0, 'max_stretch': 100.0}
[DEBUG: 2024-04-02 13:11:07 : trollimage.xrimage] Applying gamma 1.5
[DEBUG: 2024-04-02 13:11:07 : rasterio.env] Entering env context: <rasterio.env.Env object at 0x7f9b052d5fd0>
[DEBUG: 2024-04-02 13:11:07 : rasterio.env] Starting outermost env
[DEBUG: 2024-04-02 13:11:07 : rasterio.env] No GDAL environment exists
[DEBUG: 2024-04-02 13:11:07 : rasterio.env] New GDAL environment <rasterio._env.GDALEnv object at 0x7f9b04f368c0> created
[DEBUG: 2024-04-02 13:11:07 : rasterio._filepath] Installing FilePath filesystem handler plugin...
[DEBUG: 2024-04-02 13:11:07 : rasterio._env] GDAL_DATA found in environment.
[DEBUG: 2024-04-02 13:11:07 : rasterio._env] PROJ_DATA found in environment.
[DEBUG: 2024-04-02 13:11:07 : rasterio._env] Started GDALEnv: self=<rasterio._env.GDALEnv object at 0x7f9b04f368c0>.
[DEBUG: 2024-04-02 13:11:07 : rasterio.env] Entered env context: <rasterio.env.Env object at 0x7f9b052d5fd0>
[DEBUG: 2024-04-02 13:11:07 : rasterio._io] Path: _UnparsedPath(path='/tmp/test/MTG-I1-fci-202309281230-test-vis_06-test.tif'), mode: w, driver: GTiff
[DEBUG: 2024-04-02 13:11:07 : rasterio._base] Nodata success: 1, Nodata value: 0.000000
[DEBUG: 2024-04-02 13:11:07 : rasterio.env] Exiting env context: <rasterio.env.Env object at 0x7f9b052d5fd0>
[DEBUG: 2024-04-02 13:11:07 : rasterio.env] Cleared existing <rasterio._env.GDALEnv object at 0x7f9b04f368c0> options
[DEBUG: 2024-04-02 13:11:07 : rasterio._env] Stopped GDALEnv <rasterio._env.GDALEnv object at 0x7f9b04f368c0>.
[DEBUG: 2024-04-02 13:11:07 : rasterio.env] Exiting outermost env
[DEBUG: 2024-04-02 13:11:07 : rasterio.env] Exited env context: <rasterio.env.Env object at 0x7f9b052d5fd0>
[INFO: 2024-04-02 13:11:07 : satpy.writers] Computing and writing results...

Text output of actual results or error messages including full tracebacks if applicable.

Screenshots

Gwenview shows the image sliced, with the slicing somewhat unpredictable:

In The Gimp, the image looks too dark at first. When zooming in to look at the pixels, this is revealed to be due to black lines between lines with pixels:

When I remove fill_value=0 to produce an LA image, the transparency shows through in Image Magick (shown here) and The Gimp:

When I use Image Magick to convert from .tif to .png, the black lines also appear in gwenview.

Environment Info:

• OS: openSUSE Leap 15.3
• Satpy Version: v0.47.0-71-g9952f3981
• PyResample Version: v1.28.2-6-g61fd133

Additional context

When I change the resolution in the AreaDefinition from 3000 to 9000, the resulting image is all-black.

When I change the area extent from [-10000000, -10000000, 10000000, 10000000] to [-5000000, -5000000, 5000000, 5000000], the resulting image looks correct.

The problem seems to occur if the corner pixels have no data.