satellogic / telluric

telluric is a Python library to manage vector and raster geospatial data in an interactive and easy way
MIT License
87 stars 18 forks source link

GeoRaster2 `__equal__` is not always right #79

Open arielze opened 6 years ago

arielze commented 6 years ago

There are cases when two rasters have the same information, but comparing them returns False, for example:

  1. two rasters with the same information and different CRS
  2. 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
drnextgis commented 2 years ago

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.