Chunked array leads to hole in data

[saschahofmann]

I encountered an issue where a certain combination of chunksize and final resolution is creating a different output than loaded data or a different chunking schema. I struggled to create minimal example but together with the attached netcdf the problem should be reproducible by running this:

from odc.geo.geobox import GeoBox
from odc.geo.xr import xr_reproject
from pyproj import CRS
from matplotlib import pyplot as plt

CORDEX_GRID = GeoBox.from_bbox(
    (-44.6, 21.9, 65, 72.6), CRS.from_epsg(4326), resolution=0.1
)

def attrs_to_crs(attrs: dict):
    """Convert GRIB attrs to proj4 CRS."""
    latitude_of_projection_origin = attrs["GRIB_LaDInDegrees"]
    standard_parallel = attrs["GRIB_LaDInDegrees"]
    longitude_of_central_meridian = attrs["GRIB_LoVInDegrees"]
    return CRS.from_proj4(
        f"+proj=lcc +lat_1={standard_parallel} +lat_0={latitude_of_projection_origin} +lon_0={longitude_of_central_meridian} +x_0=0 +y_0=0"
    )

def add_x_y(da: xr.DataArray):
    """Adds x-y coordinates in meters and centered at 0,0"""
    dx = float(da.attrs["GRIB_DxInMetres"])
    dy = float(da.attrs["GRIB_DyInMetres"])
    return da.assign_coords(
        x=(da.x - da.x.max() / 2) * dx, y=(da.y - da.y.max() / 2) * dy
    )

def regrid_cerra(da: xr.DataArray):
    da = add_x_y(da.assign_coords(longitude=((da.longitude + 180) % 360) - 180))
    da = da.odc.assign_crs(attrs_to_crs(da.attrs))

    return xr_reproject(
        da, CORDEX_GRID, resampling='bilinear', dst_nodata=None
    )

da = xr.open_dataset('pr.nc').pr
fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))

regrid_cerra(da).plot(ax=ax1)
regrid_cerra(da.chunk(x=256, y=256)).plot(ax=ax2)

The output should look like this:

You can see that the right plot contains a small hole of NaNs. Other combinations of chunk sizes dont produce this and the in-memory example on the left also doesnt have this artifact.

[saschahofmann]

Of course, forgot to attach the file. Here you go.

pr.nc.zip

[saschahofmann]

Using different chunking removes the hole of NaNs but still introduces differences. E.g. this plot shows the differences between an in-memory version and a dask array that only has a sinlge chunk. While most of the pixels are the same there are a few "smidges" that can show up to 2°C difference.

(regrid_cerra(da)-regrid_cerra(da.chunk(x=-1, y=-1))).plot()

[Kirill888]

Thanks for the detailed report @saschahofmann

I can confirm that hole in data is a real issue and is caused by odc.geo.geobox.GeoboxTiles.{grid_intersect|tiles} methods. Error is this: output grid cell 0,2 excludes source grid cell 1,3, but should include it as they do intersect actually, but only slightly.

I think code here:

odc-geo/odc/geo/geobox.py

Lines 1584 to 1586 in 0e02b0b

	for idx in xy_chunks_with_data:
	geobox = self[idx]
	deps[idx] = list(src.tiles(geobox.extent))

should use geobox.footprint(src.base.crs, 1) instead of geobox.extent, but this is costlier in terms of runtime and also adds other issues like dateline discontinuities. What we really need is a robust yet efficient method that answers the following: "given two polygons|bounding boxes in different projections determine if they share any points in common".

Pixel differences due to chunk sizes

That's totally expected, tiny changes in coordinate representation will cause those, especially when using bilinear, but even with nearest you'll get occasional rounding into a different pixel. The cause for it is deep inside GDAL, you can probably trigger slightly different outputs by changing internal chunking for multi threaded warp operation. Both versions are valid, but experienced different floating point precision loss for coordinates, leading to output pixel differences eventually.