gradient_search comb effect

[twak]

hi, great project, love using it - it's saved me so much time!

I'm getting a 'comb' effect using the gradient resampler. I understand this is a prototype product, but it works excellently otherwise. I haven't seen this reported elsewhere, so thought I'd report.

resampling amsr2 l1r from raw h5 to a dense grid (eg from native 5000 to 80m), I get this result

if we zoom in and pull the colour range in to the view we see horizontal combs interpolating between the samples.

amsr = 'GW1AM2_201712190356_220D_L1SGRTBR_2220220', 'GW1AM2_201712190445_220A_L1SGRTBR_2220220'
roi = 324275, -1893016 : 340542, -1881592 (epsg:3031)

[djhoese]

Thanks for filing. Could you provide the exact code you ran? Also what python packages and versions are in your environment? Are you customizing any dask settings or anything like that?

[twak]

Code follows, note custom loader for satpy (amsr l1r stuff) and slooooow download of h5 from jaxa. One the tif is cooked, open (eg in qgis) and set brightness to 30-40. Look in the top-centre.

import os
import requests
from satpy.scene import Scene
import pyproj
from pyresample.geometry import AreaDefinition
from pathlib import Path


if __name__ == "__main__":

    amsr_d = 'GW1AM2_201712190356_220D_L1SGRTBR_2220220'
    amsr_f = f'{amsr_d}.h5'
    min_x, min_y, max_x, max_y = -2233691.3868798204, -247234.76145501054, -1634713.327964282, 340288.66793063114
    res = 80 # meters per pixel

    if not os.path.exists(amsr_f):

        # download from jaxa - 50Mb / 6 minutes from Europe.
        url = f'https://gportal.jaxa.jp/download/standard/GCOM-W/GCOM-W.AMSR2/L1R/2/2017/12/{amsr_f}'

        response = requests.get(url)
        print("done!")

        with open(amsr_f, 'wb') as file:
            file.write(response.content)
        print(f"downloaded {amsr_d}")

    area_def = AreaDefinition(
        area_id=f"...",
        description=f"to match footprint",
        proj_id=f'amsr_chip_proj',
        projection='EPSG:3031',
        width=(max_x-min_x)//res, height=(max_x-min_x)//res,
        area_extent=[min_x, min_y, max_x, max_y]
    )

    dss = ['btemp_89.0av']

    # custom loader (but not code) - see https://github.com/pytroll/satpy/pull/3241
    scene = Scene(reader="amsr2_l1r", filenames=[Path(amsr_f)])
    scene.load(dss)
    resampled_scene = scene.resample(area_def, datasets=dss, resampler="gradient_search", radius_of_influence=1000000)
    resampled_scene.save_datasets(
        filename="gradient.tif",
        writer='geotiff',
        fill_value=None,
    )

[djhoese]

I assume if you use resampler="nearest" in the Scene.resample call you don't see these types of patterns?

[twak]

nearest won't run on this machine at that resolution. But the patches I've tried have been okay:

(area min_y, min_x, max_y, max_x = 292480, -1910837, 343434, -1860611)

However, bilinear does this:

I've also noticed missing pixels with gradient_search at these high resolutions.

My immediate solution is to resample twice (with gradient) first from swath 5000 to raster 2000, then a second time to 80. But this is messy and costs accuracy etc...

[mraspaud]

I’m guessing this is an issue with precision in the grid position with the gradient resampling with such a jump in resolution.
Could you try to reproduce this error with synthetic data somehow, so that I can have a look at it?

[twak]

ah, let me be more helpful - i ported the repo to L1B (so satpy supports it out of the box). Download data and unzip from here. Scroll down to Sample Data -> Level-1B.

import os
import requests
from satpy.scene import Scene
import pyproj
from pyresample.geometry import AreaDefinition
from pathlib import Path

# repro for issue: https://github.com/pytroll/pyresample/issues/682#issuecomment-3333309960


# data here: https://www.eorc.jaxa.jp/AMSR/resources/index_en.html#format
# scroll down to Sample Data -> Level-1B
if __name__ == "__main__":

    amsr_d = 'GW1AM2_201301010158_201D_L1SGBTBR_2220220'
    amsr_f = f'{amsr_d}.h5'

    min_y, min_x, max_y, max_x = 58.211462, 4.660823, 60.299158, 6.740594
    res = 8000 


    dss = ['btemp_89.0av']
    scene = Scene(reader="amsr2_l1b", filenames=[Path(amsr_f)])
    scene.load(dss)

    area_def = AreaDefinition(
        area_id=f"...",
        description=f"to match footprint",
        proj_id=f'amsr_chip_proj',
        projection='EPSG:4326',
        width=(max_x-min_x) * res, height=(max_x-min_x) * res,
        area_extent=[min_x, min_y, max_x, max_y]
    )

    resampled_scene = scene.resample(area_def, datasets=dss, resampler="gradient_search", radius_of_influence=1000000)
    resampled_scene.save_datasets(
        filename="gradient.tif",
        writer='geotiff',
        fill_value=None,
    )

[mraspaud]

@twak thanks for the minimal example and the data. I will try to have a look at this next week.