telluric
telluric copied to clipboard
GeoRaster2 `__equal__` is not always right
There are cases when two rasters have the same information, but comparing them returns False
, for example:
- two rasters with the same information and different CRS
- two rasters with opposite affine scale and transformed image data
Here is a code example:
import telluric as tl
import numpy as np
from affine import Affine
from telluric.constants import WEB_MERCATOR_CRS
from shapely.geometry import Point, Polygon
array = np.arange(0, 20).reshape(1, 4, 5)
array2 = np.arange(19, -1, -1).reshape(1, 4, 5)
array2.sort()
image1 = np.ma.array(array, mask=False)
image2 = np.ma.array(array2, mask=False)
aff2 = Affine.translation(0, -8) * Affine.scale(2, 2)
aff = Affine.scale(2, -2)
r1 = tl.GeoRaster2(image=image1, affine=aff, crs=WEB_MERCATOR_CRS)
r2 = tl.GeoRaster2(image=image2, affine=aff2, crs=WEB_MERCATOR_CRS)
r1 == r2 # should be true in my opinion
r1.to_png() == r2.to_png() # is True
# comparing different points of the rasters shows they are the same
np.array_equal(r1.get(tl.GeoVector(Point(2, -2), crs=WEB_MERCATOR_CRS)),
r2.get(tl.GeoVector(Point(2, -2), crs=WEB_MERCATOR_CRS))) # is True
np.array_equal(r1.get(tl.GeoVector(Point(0, -0), crs=WEB_MERCATOR_CRS)),
r2.get(tl.GeoVector(Point(0, -0), crs=WEB_MERCATOR_CRS))) # is True
np.array_equal(r1.get(tl.GeoVector(Point(9, -2), crs=WEB_MERCATOR_CRS)),
r2.get(tl.GeoVector(Point(9, -2), crs=WEB_MERCATOR_CRS))) # is True
np.array_equal(r1.get(tl.GeoVector(Point(9, -7), crs=WEB_MERCATOR_CRS)),
r2.get(tl.GeoVector(Point(9, -7), crs=WEB_MERCATOR_CRS))) # is True
roi = tl.GeoVector(Polygon.from_bounds(0, 0, 3, -3), crs=WEB_MERCATOR_CRS)
r1c = r1.crop(roi)
r2c = r2.crop(roi)
r1c == r2c # should be true
r2c.to_png() == r1c.to_png() # is True
The fact that r2c.to_png() == r1c.to_png()
is just a coincidence. The default value of in_range
parameter of to_png
method is dtype
so it is equivalent to:
>> r1c.to_png(in_range='dtype') == r2c.to_png(in_range='dtype')
True
to_png
method calls internally astype
, but if we look closer at the images we are getting with astype
:
>>> r1c.astype(np.uint8, in_range='dtype').image
masked_array(
data=[[[127, 127],
[127, 127]]],
mask=False,
fill_value=63,
dtype=uint8)
>>>
>>> r2c.astype(np.uint8, in_range='dtype').image
masked_array(
data=[[[127, 127],
[127, 127]]],
mask=False,
fill_value=63,
dtype=uint8)
we can see that using astype
with in_range='dtype'
converts all pixel values to the same value (127).
The correct way to compare two images is (which as expected returns False
):
>>> r2c.to_png(in_range='image') == r1c.to_png(in_range='image')
False
It turn out that we still lack of the way to compare two georasters with flipped images.