Issue #761 better diff for apex reference check (replaced convex hull with bbox, _compare_xarray_dataarray_xy now only adds to the original xarray implementation)

Author: dsamaey

openeo/testing/results.py

@@ -11,18 +11,19 @@
 import numpy as np
 import xarray
 import xarray.testing
-from scipy.spatial import ConvexHull
 from xarray import DataArray
 
 from openeo.rest.job import DEFAULT_JOB_RESULTS_FILENAME, BatchJob, JobResults
 from openeo.util import repr_truncate
 
 _log = logging.getLogger(__name__)
 
-
 _DEFAULT_RTOL = 1e-6
 _DEFAULT_ATOL = 1e-6
 
+# https://paulbourke.net/dataformats/asciiart
+DEFAULT_GRAYSCALE_70_CHARACTERS = "$@B%8&WM#*oahkbdpqwmZO0QLCJUYXzcvunxrjft/\|()1{}[]?-_+~<>i!lI;:,\"^`'. "[::-1]
+DEFAULT_GRAYSCALE_10_CHARACTERS = " .:-=+*#%@"
 
 def _load_xarray_netcdf(path: Union[str, Path], **kwargs) -> xarray.Dataset:
     """
@@ -91,29 +92,33 @@ def _as_xarray_dataarray(data: Union[str, Path, xarray.DataArray]) -> xarray.Dat
     return data
 
 
-def ascii_art(diff_data: DataArray, *, max_width: int = 60, y_vs_x_aspect_ratio=2.5) -> str:
+def _ascii_art(
+    diff_data: DataArray,
+    *,
+    max_width: int = 60,
+    y_vs_x_aspect_ratio=2.5,
+    grayscale_characters: str = DEFAULT_GRAYSCALE_70_CHARACTERS,
+) -> str:
+    max_grayscale_idx = len(grayscale_characters) - 1
     x_scale: int = max(1, int(diff_data.sizes["x"] / max_width))
-    y_scale: int = max(1, int(diff_data.sizes["x"] / (max_width * y_vs_x_aspect_ratio)))
+    y_scale: int = max(1, int(diff_data.sizes["y"] / (max_width / y_vs_x_aspect_ratio)))
     data_max = diff_data.max().item()
     if data_max == 0:
         data_max = 1
-    grayscale_characters = "$@B%8&WM#*oahkbdpqwmZO0QLCJUYXzcvunxrjft/\|()1{}[]?-_+~<>i!lI;:,\"^`'. "
     coarsened = diff_data.coarsen(dim={"x": x_scale, "y": y_scale}, boundary="pad").mean()
     coarsened = coarsened.transpose("y", "x", ...)
     top = "┌" + "─" * coarsened.sizes["x"] + "┐\n"
     bottom = "\n└" + "─" * coarsened.sizes["x"] + "┘"
 
-    def pixelChar(v) -> str:
-        if np.isnan(v):
-            return " "
-        i = int(v * 69 / data_max)
+    def _pixel_char(v) -> str:
+        i = 0 if np.isnan(v) else int(v * max_grayscale_idx / data_max)
         if v > 0 and i == 0:
-            i = 1
+            i = 1  # don't show a blank for a difference above the threshold
         else:
-            i = min(69, i)
-        return grayscale_characters[69 - i]
+            i = min(max_grayscale_idx, i)
+        return grayscale_characters[i]
 
-    return top + "\n".join(["│" + "".join([pixelChar(v) for v in row]) + "│" for row in coarsened]) + bottom
+    return top + "\n".join(["│" + "".join([_pixel_char(v) for v in row]) + "│" for row in coarsened]) + bottom
 
 
 def _compare_xarray_dataarray_xy(
@@ -122,96 +127,60 @@ def _compare_xarray_dataarray_xy(
     *,
     rtol: float = _DEFAULT_RTOL,
     atol: float = _DEFAULT_ATOL,
+    name: str = None,
 ) -> List[str]:
     """
-    Compare two xarray DataArrays with tolerance and report mismatch issues (as strings)
-
-    Checks that are done (with tolerance):
-    - (optional) Check fraction of mismatching pixels (difference exceeding some tolerance).
-      If fraction is below a given threshold, ignore these mismatches in subsequent comparisons.
-      If fraction is above the threshold, report this issue.
-    - Compare actual and expected data with `xarray.testing.assert_allclose` and specified tolerances.
-
+    Additional compare for two compatible spatial xarray DataArrays with tolerance (rtol, atol)
     :return: list of issues (empty if no issues)
     """
-    # TODO: make this a public function?
-    # TODO: option for nodata fill value?
-    # TODO: option to include data type check?
-    # TODO: option to cast to some data type (or even rescale) before comparison?
-    # TODO: also compare attributes of the DataArray?
-    actual = _as_xarray_dataarray(actual)
-    expected = _as_xarray_dataarray(expected)
     issues = []
-
-    if actual.dims != expected.dims:
-        issues.append(f"Dimension mismatch: {actual.dims} != {expected.dims}")
-    for dim in sorted(set(expected.dims).intersection(actual.dims)):
-        acs = actual.coords[dim].values
-        ecs = expected.coords[dim].values
-        if not (acs.shape == ecs.shape and (acs == ecs).all()):
-            issues.append(f"Coordinates mismatch for dimension {dim!r}: {acs} != {ecs}")
-    if actual.shape != expected.shape:
-        issues.append(f"Shape mismatch: {actual.shape} != {expected.shape}")
-
-    if not issues:
-        threshold = abs(expected * rtol) + atol
-        diff_exact = abs(expected - actual)
-        diff_mask = diff_exact > threshold
-        diff_lenient = diff_exact.where(diff_mask)
-
-        non_x_y_dims = list(set(expected.dims) - {"x", "y"})
-        value_mapping = dict(map(lambda d: (d, expected[d].data), non_x_y_dims))
-        shape = tuple([len(value_mapping[x]) for x in non_x_y_dims])
-
-        for shape_index, v in np.ndenumerate(np.ndarray(shape)):
-            indexers = {}
-            for index, value_index in enumerate(shape_index):
-                indexers[non_x_y_dims[index]] = value_mapping[non_x_y_dims[index]][value_index]
-            diff_data = diff_lenient.sel(indexers=indexers)
-            total_pixel_count = expected.sel(indexers).count().item()
-            diff_pixel_count = diff_data.count().item()
-
-            if diff_pixel_count > 0:
-                diff_pixel_percentage = round(diff_pixel_count * 100 / total_pixel_count, 1)
-                diff_mean = round(diff_data.mean().item(), 2)
-                diff_var = round(diff_data.var().item(), 2)
-
-                key = ",".join([f"{k} {str(v1)}" for k, v1 in indexers.items()])
-                issues.append(
-                    f"{key}: value difference exceeds tolerance (rtol {rtol}, atol {atol}), min:{diff_data.min().data}, max: {diff_data.max().data}, mean: {diff_mean}, var: {diff_var}"
-                )
-
-                coord_grid = np.meshgrid(diff_data.coords["x"], diff_data.coords["y"])
-
-                mask = diff_data.notnull()
-                if mask.dims[0] != "y":
-                    mask = mask.transpose()
-                c1 = coord_grid[0][mask]
-                c2 = coord_grid[1][mask]
-                coordinates = np.dstack((c1, c2)).reshape(-1, 2)
-
-                art = ascii_art(diff_data)
-                print(f"Difference ascii art for {key}")
-                print(art)
-
-                if len(coordinates) > 2:
-                    hull = ConvexHull(coordinates)
-                    area = hull.volume
-
-                    x_m = diff_data.coords["x"][0].data
-                    x_M = diff_data.coords["x"][-1].data
-                    y_m = diff_data.coords["y"][0].data
-                    y_M = diff_data.coords["y"][-1].data
-
-                    total_area = abs((y_M - y_m) * (x_M - x_m))
-                    area_percentage = round(area * 100 / total_area, 1)
-                    issues.append(
-                        f"{key}: differing pixels: {diff_pixel_count}/{total_pixel_count} ({diff_pixel_percentage}%), spread over {area_percentage}% of the area"
-                    )
-                else:
-                    issues.append(
-                        f"{key}: differing pixels: {diff_pixel_count}/{total_pixel_count} ({diff_pixel_percentage}%)"
-                    )
+    threshold = abs(expected * rtol) + atol
+    diff_exact = abs(expected - actual)
+    diff_mask = diff_exact > threshold
+    diff_lenient = diff_exact.where(diff_mask)
+
+    non_x_y_dims = list(set(expected.dims) - {"x", "y"})
+    value_mapping = dict(map(lambda d: (d, expected[d].data), non_x_y_dims))
+    shape = tuple([len(value_mapping[x]) for x in non_x_y_dims])
+
+    for shape_index, v in np.ndenumerate(np.ndarray(shape)):
+        indexers = {}
+        for index, value_index in enumerate(shape_index):
+            indexers[non_x_y_dims[index]] = value_mapping[non_x_y_dims[index]][value_index]
+        diff_data = diff_lenient.sel(indexers=indexers)
+        total_pixel_count = expected.sel(indexers).count().item()
+        diff_pixel_count = diff_data.count().item()
+
+        if diff_pixel_count > 0:
+            diff_pixel_percentage = round(diff_pixel_count * 100 / total_pixel_count, 1)
+            diff_mean = round(diff_data.mean().item(), 2)
+            diff_var = round(diff_data.var().item(), 2)
+
+            key = name + ": " if name else ""
+            key += ",".join([f"{k} {str(v1)}" for k, v1 in indexers.items()])
+            issues.append(
+                f"{key}: value difference exceeds tolerance (rtol {rtol}, atol {atol}), min:{diff_data.min().data}, max: {diff_data.max().data}, mean: {diff_mean}, var: {diff_var}"
+            )
+
+            _log.warning(f"Difference (ascii art) for {key}:\n{_ascii_art(diff_data)}")
+
+            coord_grid = np.meshgrid(diff_data.coords["x"], diff_data.coords["y"])
+            mask = diff_data.notnull()
+            if mask.dims[0] != "y":
+                mask = mask.transpose()
+            x_coords = coord_grid[0][mask]
+            y_coords = coord_grid[1][mask]
+
+            diff_bbox = ((x_coords.min().item(), y_coords.min().item()), (x_coords.max().item(), y_coords.max().item()))
+            diff_area = (x_coords.max() - x_coords.min()) * (y_coords.max() - y_coords.min())
+            total_area = abs(
+                (diff_data.coords["y"][-1].data - diff_data.coords["y"][0].data)
+                * (diff_data.coords["x"][-1].data - diff_data.coords["x"][0].data)
+            )
+            area_percentage = round(diff_area * 100 / total_area, 1)
+            issues.append(
+                f"{key}: differing pixels: {diff_pixel_count}/{total_pixel_count} ({diff_pixel_percentage}%), bbox {diff_bbox} - {area_percentage}% of the area"
+            )
     return issues
 
 
@@ -221,6 +190,7 @@ def _compare_xarray_dataarray(
     *,
     rtol: float = _DEFAULT_RTOL,
     atol: float = _DEFAULT_ATOL,
+    name: str = None,
 ) -> List[str]:
     """
     Compare two xarray DataArrays with tolerance and report mismatch issues (as strings)
@@ -243,7 +213,7 @@ def _compare_xarray_dataarray(
     issues = []
 
     # `xarray.testing.assert_allclose` currently does not always
-    # provides detailed information about shape/dimension mismatches
+    # provide detailed information about shape/dimension mismatches
     # so we enrich the issue listing with some more details
     if actual.dims != expected.dims:
         issues.append(f"Dimension mismatch: {actual.dims} != {expected.dims}")
@@ -254,17 +224,14 @@ def _compare_xarray_dataarray(
             issues.append(f"Coordinates mismatch for dimension {dim!r}: {acs} != {ecs}")
     if actual.shape != expected.shape:
         issues.append(f"Shape mismatch: {actual.shape} != {expected.shape}")
-
-    if not issues:
-        if {"x", "y"} <= set(expected.dims):
-            issues = _compare_xarray_dataarray_xy(actual=actual, expected=expected, rtol=rtol, atol=atol)
-        else:
-            try:
-                xarray.testing.assert_allclose(a=actual, b=expected, rtol=rtol, atol=atol)
-            except AssertionError as e:
-                # TODO: message of `assert_allclose` is typically multiline, split it again or make it one line?
-                issues.append(str(e).strip())
-
+    compatible = len(issues) == 0
+    try:
+        xarray.testing.assert_allclose(a=actual, b=expected, rtol=rtol, atol=atol)
+    except AssertionError as e:
+        # TODO: message of `assert_allclose` is typically multiline, split it again or make it one line?
+        issues.append(str(e).strip())
+    if compatible and {"x", "y"} <= set(expected.dims):
+        issues.extend(_compare_xarray_dataarray_xy(actual=actual, expected=expected, rtol=rtol, atol=atol, name=name))
     return issues
 
 
@@ -293,32 +260,6 @@ def assert_xarray_dataarray_allclose(
     if issues:
         raise AssertionError("\n".join(issues))
 
-
-def assert_xarray_dataarray_allclose_xy(
-    actual: Union[xarray.DataArray, str, Path],
-    expected: Union[xarray.DataArray, str, Path],
-    *,
-    rtol: float = _DEFAULT_RTOL,
-    atol: float = _DEFAULT_ATOL,
-):
-    """
-    Assert that two Xarray ``DataArray`` instances are equal (with tolerance).
-
-    :param actual: actual data, provided as Xarray DataArray object or path to NetCDF/GeoTIFF file.
-    :param expected: expected or reference data, provided as Xarray DataArray object or path to NetCDF/GeoTIFF file.
-    :param rtol: relative tolerance
-    :param atol: absolute tolerance
-    :raises AssertionError: if not equal within the given tolerance
-
-    .. versionadded:: 0.31.0
-
-    .. warning::
-        This function is experimental and subject to change.
-    """
-    issues = _compare_xarray_dataarray_xy(actual=actual, expected=expected, rtol=rtol, atol=atol)
-    if issues:
-        raise AssertionError("\n".join(issues))
-
 def _compare_xarray_datasets(
     actual: Union[xarray.Dataset, str, Path],
     expected: Union[xarray.Dataset, str, Path],
@@ -336,15 +277,14 @@ def _compare_xarray_datasets(
     expected = _as_xarray_dataset(expected)
 
     all_issues = []
-    # TODO: just leverage DataSet support in xarray.testing.assert_allclose for all this?
     actual_vars = set(actual.data_vars)
     expected_vars = set(expected.data_vars)
     _log.debug(f"_compare_xarray_datasets: actual_vars={actual_vars!r} expected_vars={expected_vars!r}")
     if actual_vars != expected_vars:
         all_issues.append(f"Xarray DataSet variables mismatch: {actual_vars} != {expected_vars}")
     for var in expected_vars.intersection(actual_vars):
         _log.debug(f"_compare_xarray_datasets: comparing variable {var!r}")
-        issues = _compare_xarray_dataarray(actual[var], expected[var], rtol=rtol, atol=atol)
+        issues = _compare_xarray_dataarray(actual[var], expected[var], rtol=rtol, atol=atol, name=var)
         if issues:
             all_issues.append(f"Issues for variable {var!r}:")
             all_issues.extend(issues)
@@ -406,10 +346,7 @@ def assert_xarray_allclose(
     if isinstance(actual, xarray.Dataset) and isinstance(expected, xarray.Dataset):
         assert_xarray_dataset_allclose(actual, expected, rtol=rtol, atol=atol)
     elif isinstance(actual, xarray.DataArray) and isinstance(expected, xarray.DataArray):
-        if {"x", "y"}.issubset(expected.dims):
-            assert_xarray_dataarray_allclose_xy(actual, expected, rtol=rtol, atol=atol)
-        else:
-            assert_xarray_dataarray_allclose(actual, expected, rtol=rtol, atol=atol)
+        assert_xarray_dataarray_allclose(actual, expected, rtol=rtol, atol=atol)
     else:
         raise ValueError(f"Unsupported types: {type(actual)} and {type(expected)}")
 

setup.py

@@ -82,7 +82,6 @@
         "deprecated>=1.2.12",
         'oschmod>=0.3.12; sys_platform == "win32"',
         "importlib_resources; python_version<'3.9'",
-        "scipy",  # for Convex Hull algorithm
     ],
     extras_require={
         "tests": tests_require,