We have a lot of steps that extract the NumPy array/transform from a Raster/Xarray dataset, and then re-build the object in the fit function to be able to apply the Raster API consistently. This takes extra memory and computing time.

I see 2 solutions to remedy this issue:

1. Switch to using Raster and Xarray objects directly in the fit_rst, never going transforming to a NumPy array (just passing the one contained),
2. Make the Raster API available as separate functions from the class.

We (wisely) decided to not do 1. initially because we did not know if the Raster class was going to be a success. But given that we'll be able to use an Xarray object the same way now, that should not be a problem anymore.
It would still be valuable to have 2. available for some instances.

Minimal changes would be needed for 2., but would make the code harder to read in the Coreg class (calling non-public GeoUtils functions with volatile metadata).