Pyresample with H5NetCDF and earthaccess

import earthaccess
import xarray as xr
from pyresample.area_config import create_area_def
from pyresample.gradient import block_nn_interpolator, gradient_resampler_indices_block
from pyresample.resampler import resample_blocks

def warp_resample():
    bucket = "podaac-ops-cumulus-protected"
    input_uri = "MUR-JPL-L4-GLOB-v4.1/20020601090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc"
    variable = "analysed_sst"
    src = f"s3://{bucket}/{input_uri}"
    dstSRS = "EPSG:3857"
    srcSRS = "EPSG:4326"
    width = height = 256
    te = [
        -20037508.342789244,
        -20037508.342789244,
        20037508.342789244,
        20037508.342789244,
    ]
    earthaccess.login()
    fs = earthaccess.get_s3fs_session(daac="PODAAC")
    fsspec_caching = {
        "cache_type": "none",
    }
    with fs.open(src, **fsspec_caching) as f:
        da = xr.open_dataset(f, engine="h5netcdf")[variable]
        da = da.chunk({"time": -1, "lat": 4000, "lon": 4000})
        target_area_def = create_area_def(
            area_id=1,
            projection=dstSRS,
            shape=(height, width),
            area_extent=te,
        )
        source_area_def = create_area_def(
            area_id=2,
            projection=srcSRS,
            shape=(da.sizes["lat"], da.sizes["lon"]),
            area_extent=[-179.995, 89.995, 180.005, -89.995],
        )
        indices_xy = resample_blocks(
            gradient_resampler_indices_block,
            source_area_def,
            [],
            target_area_def,
            chunk_size=(1, height, width),
            dtype=float,
        )
        resampled = resample_blocks(
            block_nn_interpolator,
            source_area_def,
            [da.data],
            target_area_def,
            dst_arrays=[indices_xy],
            chunk_size=(1, height, width),
            dtype=da.dtype,
        )
        return resampled.compute()

if __name__ == "__main__":
    warp_resample()