Discovered that xarray.stack loads the data into memory (https://github.com/pydata/xarray/issues/4113). aggregate spatial can work without flattening the data first, so this reorganises it to avoid .stack.
Also I've added a fix to make sure the mask always has the correct crs set.
Discovered that xarray.stack loads the data into memory (https://github.com/pydata/xarray/issues/4113). aggregate spatial can work without flattening the data first, so this reorganises it to avoid .stack. Also I've added a fix to make sure the mask always has the correct crs set.