Asymmetry of arithmetic operations when broadcasting between 2D and 3D arrays

[rhugonnet]

This actually works now with #302 (the numpy interface is not taking over anymore):

test = np.ones(r.shape) + r
type(test)
Out[15]: geoutils.georaster.raster.Raster
test = r + np.ones(r.data.shape)
type(test)
Out[17]: geoutils.georaster.raster.Raster

However, as you can notice, I have to change the shape to .data.shape on one side, and not another, otherwise:

ValueError: Both rasters must have the same shape.

I'll leave this issue to be closed by #302 and open a new one on this.

Originally posted by @rhugonnet in #256 (comment)