all_touched=True unexpected behavior

[ivanhigueram]

Hello all!

I am using the geocube.api.core.make_cube function as a way to build a template array for fast data extraction in my pipeline. Basically, I will use this template array to store the row-index of a pd.GeoDataFrame that overlap with each of the grids of my ds (a xr.DataArray) My prior with the all_touched rasterization process is that when True, all pixels that overlap with the raster will be included, and if False, only pixels whose centroid overlaps with the geometry will be included.

In my application, I expect that if all_touched=True, all row-indices will appear in my template, as even small geometries will overlap with any pixel, but I found the opposite.

The testing suite does not have a good nc file, so I created a random one with the soil data:

# Testing

gdf = gpd.read_file("./soil_data_flat.geojson")
gdf.reset_index(drop=False, inplace=True)

# Generate random array
# Get the bounds of the GeoDataFrame
bounds = gdf.total_bounds
minx, miny, maxx, maxy = bounds

# Generate a grid of longitude and latitude values
lon = np.linspace(minx, maxx, 100)
lat = np.linspace(miny, maxy, 100)

# Generate random data for the grid
data = np.random.rand(len(lat), len(lon))

# Create the DataArray
da = xr.DataArray(
    data,
    coords=[('lat', lat), ('lon', lon)],
    name='random_data'
)
da.rio.set_spatial_dims('lon', 'lat', inplace=True)
da.rio.set_crs(gdf.crs, inplace=True)

Now, with all_touched=False:

from geocube.api.core import make_geocube
from functools import partial
from geocube.rasterize import rasterize_image
from rasterio.enums import MergeAlg 

template = make_geocube(
    vector_data=gdf[["index", "geometry"]],
     measurements=["index"],
      like=da,
      rasterize_function=partial(rasterize_image, all_touched=False),
)

np.unique(template.index.values).shape
> (8, )

and with all_touched=True:

template = make_geocube(
    vector_data=gdf[["index", "geometry"]],
     measurements=["index"],
      like=da,
      rasterize_function=partial(rasterize_image, all_touched=True),
)

np.unique(template.index.values).shape
> (7, )

Is this behavior expected? I am misusing the make_geocube function?

[snowman2]

If you look at the values:

{
    "original": numpy.unique(original_raster.index).tolist(),
    "all_touched": numpy.unique(all_touched_raster.index).tolist(),
}

all_touched=True appears to fill in the nan values at the edges:

{'original': [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, nan],
 'all_touched': [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]}

I created a raster diff layer on top of the origin input polygons to illustrate the differences between the two and it appears to work as expected.

Difference mask:

different_mask = all_touched_raster.index != original_raster.index

Raster all_touched=False where it is different from all_touched=True overlayed on the original shapefile:

original_raster.index.where(different_mask)

In this scenario, the differing cells use the values within the shapes they touch and don't extend past the boundaries.

Raster all_touched=True where it is different from all_touched=False overlayed on the original shapefile:

all_touched_raster.index.where(different_mask)

In this scenario, the values in the cells often are pulled from the shapes next to the cells as well as values on the boundaries of the raster.