Add custom exception when Raster (sub)classes are supplied to Coreg classes

[erikmannerfelt]

I quite often accidentally do:

coreg = xdem.coreg.ICP()

reference_dem = xdem.DEM(...)

coreg.fit(reference_dem, ...) 

note that this will not work, as the reference_dem argument is supposed to be a numpy array (so the function call should take reference_dem.data). When the above is called, a not-so-helpful error message comes up:

Traceback (most recent call last):
  File "/home/erik/Projects/ETH/xdem/hello.py", line 36, in <module>
    dem_coreg = piecewise2.apply(dem_1990, transform=dem_1990.transform)
  File "/home/erik/Projects/ETH/xdem/xdem/coreg.py", line 534, in apply
    dem_array, dem_mask = xdem.spatial_tools.get_array_and_mask(dem)
  File "/home/erik/Projects/ETH/xdem/xdem/spatial_tools.py", line 47, in get_array_and_mask
    invalid_mask = get_mask(array)
  File "/home/erik/Projects/ETH/xdem/xdem/spatial_tools.py", line 25, in get_mask
    mask = (array.mask | ~np.isfinite(array.data)) if isinstance(array, np.ma.masked_array) else ~np.isfinite(array)
TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

I see three solutions for this:

1. Add a custom exception for isinstance(reference_dem, Raster) to be more helpful.
2. Allow Raster (sub)classes as coreg arguments. This would reduce the need for the transform argument as that can be read automatically from reference_dem.
3. Allow functions like np.isfinite to run on a Raster as if it was called on Raster.data: np.isfinite(Raster) -> np.ndarray

Thoughs, @adehecq and @rhugonnet ?

[adehecq]

• Point 3 would be great, but it's a lot more difficult than it seems. If you have any idea on how to do this, it would be great (I guess Pandas dataframe behave this way?).
• Point 2 raises again the question whether we want to make xdem fully dependent on geoutils. I think it makes sense to at least accept numpy arrays in all cases. I guess an "hybrid" approach where one can provide a Raster object or (a numpy array and a transform) would be great. Then transform could be set as optional and would be used only for a numpy array. But this would need to be handled the same way in all xdem functions.
  Not sure it gives you the answer you're looking for...

[erikmannerfelt]

* Point 3 would be great, but it's a lot more difficult than it seems. If you have any idea on how to do this, it would be great (I guess Pandas dataframe behave this way?).

* Point 2 raises again the question whether we want to make xdem fully dependent on geoutils. I think it makes sense to at least accept numpy arrays in all cases. I guess an "hybrid" approach where one can provide a Raster object or (a numpy array and a transform) would be great. Then transform could be set as optional and would be used only for a numpy array. But this would need to be handled the same way in all xdem functions.
  Not sure it gives you the answer you're looking for...

We could have a raster_to_array_and_transform function that is called whenever Rasters are allowed!

[adehecq]

We could have a raster_to_array_and_transform function that is called whenever Rasters are allowed!

Well, this seems straightforward enough to not have a function for it. But you could have a condition at the beginning of the function, to define both terms either from the Raster attributes, or from the arguments raster and transform.

[rhugonnet]

We had this in GeoImg, basically each function started with a:

if isinstance(reference_dem,GeoImg):
    reference_array=reference_dem.img
elif isinstance(reference_dem,np.ndarray):
    reference_array=reference_dem
else:
    raise ValueError()

But this rapidly becomes annoying... we get used to it and, then we need to write it at the beginning of each function.

So I agree: 3 would be ideal... but seems complex.
Otherwise for 2: I guess that, using the Coreg class level, this would be less annoying? (don't need to repeat for every Coreg subclass). But once we want to do the same outside this class, we'll face the same problem again... :/

[erikmannerfelt]

Since the latest geoutils, numpy operations now work natively on rasters! This issue might therefore be solved (as soon as it is included in a release of course). But the behaviour should optimally be tested in xdem.