Tests improvements (#636)

* increase line length to 100

* generalize tests

* test aspect

* added conftest.py

* assert numpy equal dask cupy case

* test classify

* curvature to return float32

* test curvature

* test dtype of output

* classify to return float32

* set verify_dtype to False as default

* test focal

* test hillshade

* aspect, classify, curvature: verify dtype of output array

* test multispectral

* linting

* added rtol as a param

* test pathfinding

* allocation: return dtype of float32

* test proximity

* flake8

* test slope

* test terrain

* test viewshed

* test zonal

* change as requested

* Fix cupy and rtxpy tests

Co-authored-by: Ian Thomas <[email protected]>

Author: thuydotm

benchmarks/asv.conf.json

@@ -73,7 +73,7 @@
     "matrix": {
         "pyct": [],
         // Optional for cupy support; use correct cuda version.
-        //"cupy-cuda101": [],
+        //"cupy-cuda113": [],
         // Optional for ray-tracing support, need cupy too.
         //"rtxpy": [],
         // Optional for polygonize.

xrspatial/classify.py

@@ -243,7 +243,7 @@ def reclassify(agg: xr.DataArray,
         >>> agg_da = xr.DataArray(data_da, name='agg_da')
         >>> print(agg_da)
         <xarray.DataArray 'agg_da' (dim_0: 4, dim_1: 5)>
-        dask.array<array, shape=(4, 5), dtype=float64, chunksize=(3, 3), chunktype=numpy.ndarray>
+        dask.array<array, shape=(4, 5), dtype=float32, chunksize=(3, 3), chunktype=numpy.ndarray>
         Dimensions without coordinates: dim_0, dim_1
         >>> agg_reclassify_da = reclassify(agg_da, bins=bins, new_values=new_values)  # noqa
         >>> print(agg_reclassify_da)
@@ -402,12 +402,12 @@ def quantile(agg: xr.DataArray,
 def _run_numpy_jenks_matrices(data, n_classes):
     n_data = data.shape[0]
     lower_class_limits = np.zeros(
-        (n_data + 1, n_classes + 1), dtype=np.float64
+        (n_data + 1, n_classes + 1), dtype=np.float32
     )
     lower_class_limits[1, 1:n_classes + 1] = 1.0
 
     var_combinations = np.zeros(
-        (n_data + 1, n_classes + 1), dtype=np.float64
+        (n_data + 1, n_classes + 1), dtype=np.float32
     )
     var_combinations[2:n_data + 1, 1:n_classes + 1] = np.inf
 
@@ -423,7 +423,7 @@ def _run_numpy_jenks_matrices(data, n_classes):
             lower_class_limit = l - m
             i4 = lower_class_limit - 1
 
-            val = np.float64(data[i4])
+            val = np.float32(data[i4])
 
             # here we're estimating variance for each potential classing
             # of the data, for each potential number of classes. `w`
@@ -521,7 +521,7 @@ def _run_jenks(data, n_classes, module):
         lower_class_limits, _ = _run_cupy_jenks_matrices(data, n_classes)
 
     k = data.shape[0]
-    kclass = np.zeros(n_classes + 1, dtype=np.float64)
+    kclass = np.zeros(n_classes + 1, dtype=np.float32)
     kclass[0] = data[0]
     kclass[-1] = data[-1]
     count_num = n_classes

xrspatial/curvature.py

@@ -27,7 +27,7 @@ class cupy(object):
 
 @ngjit
 def _cpu(data, cellsize):
-    out = np.empty(data.shape, np.float64)
+    out = np.empty(data.shape, np.float32)
     out[:] = np.nan
     rows, cols = data.shape
     for y in range(1, rows - 1):
@@ -128,7 +128,7 @@ def curvature(agg: xr.DataArray,
         >>> import numpy as np
         >>> import xarray as xr
         >>> from xrspatial import curvature
-        >>> flat_data = np.zeros((5, 5), dtype=np.float64)
+        >>> flat_data = np.zeros((5, 5), dtype=np.float32)
         >>> flat_raster = xr.DataArray(flat_data, attrs={'res': (1, 1)})
         >>> flat_curv = curvature(flat_raster)
         >>> print(flat_curv)
@@ -149,20 +149,20 @@ def curvature(agg: xr.DataArray,
             [0, 0, 0, 0, 0],
             [0, 0, -1, 0, 0],
             [0, 0, 0, 0, 0],
-            [0, 0, 0, 0, 0]], dtype=np.float64)
+            [0, 0, 0, 0, 0]], dtype=np.float32)
         >>> convex_raster = xr.DataArray(
             da.from_array(convex_data, chunks=(3, 3)),
             attrs={'res': (10, 10)}, name='convex_dask_numpy_raster')
         >>> print(convex_raster)
         <xarray.DataArray 'convex_dask_numpy_raster' (dim_0: 5, dim_1: 5)>
-        dask.array<array, shape=(5, 5), dtype=float64, chunksize=(3, 3), chunktype=numpy.ndarray>
+        dask.array<array, shape=(5, 5), dtype=float32, chunksize=(3, 3), chunktype=numpy.ndarray>
         Dimensions without coordinates: dim_0, dim_1
         Attributes:
             res:      (10, 10)
         >>> convex_curv = curvature(convex_raster, name='convex_curvature')
         >>> print(convex_curv)  # return a xarray DataArray with Dask-backed array
         <xarray.DataArray 'convex_curvature' (dim_0: 5, dim_1: 5)>
-        dask.array<_trim, shape=(5, 5), dtype=float64, chunksize=(3, 3), chunktype=numpy.ndarray>
+        dask.array<_trim, shape=(5, 5), dtype=float32, chunksize=(3, 3), chunktype=numpy.ndarray>
         Dimensions without coordinates: dim_0, dim_1
         Attributes:
             res:      (10, 10)
@@ -185,7 +185,7 @@ def curvature(agg: xr.DataArray,
             [0, 0, 0, 0, 0],
             [0, 0, 1, 0, 0],
             [0, 0, 0, 0, 0],
-            [0, 0, 0, 0, 0]], dtype=np.float64)
+            [0, 0, 0, 0, 0]], dtype=np.float32)
         >>> concave_raster = xr.DataArray(
             cupy.asarray(concave_data),
             attrs={'res': (10, 10)}, name='concave_cupy_raster')

xrspatial/proximity.py

@@ -885,7 +885,7 @@ def allocation(
 
     # convert to have same type as of input @raster
     result = xr.DataArray(
-        (allocation_img).astype(raster.dtype),
+        allocation_img,
         coords=raster.coords,
         dims=raster.dims,
         attrs=raster.attrs,

xrspatial/tests/conftest.py

@@ -0,0 +1,11 @@
+import pytest
+
+import numpy as np
+
+
+@pytest.fixture
+def random_data(size, dtype):
+    rng = np.random.default_rng(2841)
+    data = rng.integers(-100, 100, size=size)
+    data = data.astype(dtype)
+    return data

xrspatial/tests/general_checks.py

@@ -2,13 +2,32 @@
 import dask.array as da
 import xarray as xr
 
+from xrspatial.utils import has_cuda
 from xrspatial.utils import ArrayTypeFunctionMapping
 
 
+def create_test_raster(data, backend='numpy', dims=['y', 'x'], attrs=None, chunks=(3, 3)):
+    raster = xr.DataArray(data, dims=dims, attrs=attrs)
+    # set coords for test raster
+    for i, dim in enumerate(dims):
+        raster[dim] = np.linspace(0, data.shape[i] - 1, data.shape[i])
+
+    if has_cuda() and 'cupy' in backend:
+        import cupy
+        raster.data = cupy.asarray(raster.data)
+
+    if 'dask' in backend:
+        raster.data = da.from_array(raster.data, chunks=chunks)
+
+    return raster
+
+
 def general_output_checks(input_agg: xr.DataArray,
                           output_agg: xr.DataArray,
                           expected_results: np.ndarray = None,
-                          verify_attrs: bool = True):
+                          verify_attrs: bool = True,
+                          verify_dtype: bool = False,
+                          rtol=1e-06):
 
     # type of output is the same as of input
     assert isinstance(output_agg.data, type(input_agg.data))
@@ -29,23 +48,64 @@ def general_output_checks(input_agg: xr.DataArray,
             )
 
     if expected_results is not None:
-        numpy_func = lambda output, expected: np.testing.assert_allclose(  # noqa: E731, E501
-            output, expected_results, equal_nan=True, rtol=1e-06
-        )
-        dask_func = lambda output, expected: np.testing.assert_allclose(  # noqa: E731, E501
-            output.compute(), expected_results, equal_nan=True, rtol=1e-06
-        )
-        cupy_func = lambda output, expected: np.testing.assert_allclose(  # noqa: E731, E501
-            output.get(), expected_results, equal_nan=True, rtol=1e-06
-        )
-        dask_cupy_func = lambda output, expected: np.testing.assert_allclose(  # noqa: E731, E501
-            output.compute().get(), expected_results,
-            equal_nan=True, rtol=1e-06
-        )
+        get_numpy_data = lambda output: output  # noqa: E731
+        get_dask_numpy_data = lambda output: output.compute()  # noqa: E731
+        get_cupy_data = lambda output: output.get()  # noqa: E731
+        get_dask_cupy_data = lambda output: output.compute().get()  # noqa: E731
+
         mapper = ArrayTypeFunctionMapping(
-            numpy_func=numpy_func,
-            dask_func=dask_func,
-            cupy_func=cupy_func,
-            dask_cupy_func=dask_cupy_func,
+            numpy_func=get_numpy_data,
+            dask_func=get_dask_numpy_data,
+            cupy_func=get_cupy_data,
+            dask_cupy_func=get_dask_cupy_data,
         )
-        mapper(output_agg)(output_agg.data, expected_results)
+        output_data = mapper(output_agg)(output_agg.data)
+        np.testing.assert_allclose(output_data, expected_results, equal_nan=True, rtol=rtol)
+
+        if verify_dtype:
+            assert output_data.dtype == expected_results.dtype
+
+
+def assert_nan_edges_effect(result_agg):
+    # nan edge effect
+    edges = [
+        result_agg.data[0, :],
+        result_agg.data[-1, :],
+        result_agg.data[:, 0],
+        result_agg.data[:, -1],
+    ]
+    for edge in edges:
+        np.testing.assert_array_equal(edge, np.nan)
+
+
+def assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, func, nan_edges=True):
+    numpy_result = func(numpy_agg)
+    if nan_edges:
+        assert_nan_edges_effect(numpy_result)
+
+    dask_result = func(dask_agg)
+    general_output_checks(dask_agg, dask_result)
+    np.testing.assert_allclose(numpy_result.data, dask_result.data.compute(), equal_nan=True)
+
+
+def assert_numpy_equals_cupy(numpy_agg, cupy_agg, func, nan_edges=True, atol=0, rtol=1e-7):
+    numpy_result = func(numpy_agg)
+    if nan_edges:
+        assert_nan_edges_effect(numpy_result)
+
+    cupy_result = func(cupy_agg)
+    general_output_checks(cupy_agg, cupy_result)
+    np.testing.assert_allclose(
+        numpy_result.data, cupy_result.data.get(), equal_nan=True, atol=atol, rtol=rtol)
+
+
+def assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, func, nan_edges=True):
+    numpy_result = func(numpy_agg)
+    if nan_edges:
+        assert_nan_edges_effect(numpy_result)
+
+    dask_cupy_result = func(dask_cupy_agg)
+    general_output_checks(dask_cupy_agg, dask_cupy_result)
+    np.testing.assert_allclose(
+        numpy_result.data, dask_cupy_result.data.compute().get(), equal_nan=True
+    )

xrspatial/tests/test_aspect.py

@@ -1,92 +1,91 @@
-import dask.array as da
 import numpy as np
 import pytest
-import xarray as xr
 
 from xrspatial import aspect
 from xrspatial.utils import doesnt_have_cuda
 
+from xrspatial.tests.general_checks import create_test_raster
+from xrspatial.tests.general_checks import assert_numpy_equals_dask_numpy
+from xrspatial.tests.general_checks import assert_numpy_equals_cupy
+from xrspatial.tests.general_checks import assert_nan_edges_effect
 from xrspatial.tests.general_checks import general_output_checks
 
 
-INPUT_DATA = np.asarray([
-    [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-    [1584.8767, 1584.8767, 1585.0546, 1585.2324, 1585.2324, 1585.2324],
-    [1585.0546, 1585.0546, 1585.2324, 1585.588, 1585.588, 1585.588],
-    [1585.2324, 1585.4102, 1585.588, 1585.588, 1585.588, 1585.588],
-    [1585.588, 1585.588, 1585.7659, 1585.7659, 1585.7659, 1585.7659],
-    [1585.7659, 1585.9437, 1585.7659, 1585.7659, 1585.7659, 1585.7659],
-    [1585.9437, 1585.9437, 1585.9437, 1585.7659, 1585.7659, 1585.7659]],
-    dtype=np.float32
-)
-
-QGIS_OUTPUT = np.asarray([
-    [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-    [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-    [330.94687, 335.55496, 320.70786, 330.94464, 0., 0.],
-    [333.43494, 333.43494, 329.03394, 341.56897, 0., 18.434948],
-    [338.9621, 338.20062, 341.56506, 0., 0., 45.],
-    [341.56506, 351.8699, 26.56505, 45., -1., 90.],
-    [351.86676, 11.306906, 45., 45., 45., 108.431015]], dtype=np.float32
-)
-
+def input_data(backend='numpy'):
+    data = np.asarray([
+        [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
+        [1584.8767, 1584.8767, 1585.0546, 1585.2324, 1585.2324, 1585.2324],
+        [1585.0546, 1585.0546, 1585.2324, 1585.588, 1585.588, 1585.588],
+        [1585.2324, 1585.4102, 1585.588, 1585.588, 1585.588, 1585.588],
+        [1585.588, 1585.588, 1585.7659, 1585.7659, 1585.7659, 1585.7659],
+        [1585.7659, 1585.9437, 1585.7659, 1585.7659, 1585.7659, 1585.7659],
+        [1585.9437, 1585.9437, 1585.9437, 1585.7659, 1585.7659, 1585.7659]],
+        dtype=np.float32
+    )
+    raster = create_test_raster(data, backend, attrs={'res': (10.0, 10.0)})
+    return raster
+
+
+@pytest.fixture
+def qgis_output():
+    result = np.array([
+        [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
+        [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
+        [330.94687, 335.55496, 320.70786, 330.94464, 0., 0.],
+        [333.43494, 333.43494, 329.03394, 341.56897, 0., 18.434948],
+        [338.9621, 338.20062, 341.56506, 0., 0., 45.],
+        [341.56506, 351.8699, 26.56505, 45., -1., 90.],
+        [351.86676, 11.306906, 45., 45., 45., 108.431015]], dtype=np.float32
+    )
+    return result
 
-def test_numpy_equals_qgis():
 
-    small_da = xr.DataArray(INPUT_DATA, attrs={'res': (10.0, 10.0)})
-    xrspatial_aspect = aspect(small_da, name='numpy_aspect')
+def test_numpy_equals_qgis(qgis_output):
+    numpy_agg = input_data()
+    xrspatial_aspect = aspect(numpy_agg, name='numpy_aspect')
 
-    general_output_checks(small_da, xrspatial_aspect)
+    general_output_checks(numpy_agg, xrspatial_aspect, verify_dtype=True)
     assert xrspatial_aspect.name == 'numpy_aspect'
 
-    # validate output values
     xrspatial_vals = xrspatial_aspect.data[1:-1, 1:-1]
-    qgis_vals = QGIS_OUTPUT[1:-1, 1:-1]
+    qgis_vals = qgis_output[1:-1, 1:-1]
     # aspect is nan if nan input
     # aspect is invalid (-1) if slope equals 0
-    # otherwise aspect are from 0 - 360
+    # otherwise aspect are from 0 to 360
     np.testing.assert_allclose(xrspatial_vals, qgis_vals, equal_nan=True)
-
     # nan edge effect
-    xrspatial_edges = [
-        xrspatial_aspect.data[0, :],
-        xrspatial_aspect.data[-1, :],
-        xrspatial_aspect.data[:, 0],
-        xrspatial_aspect.data[:, -1],
-    ]
-    for edge in xrspatial_edges:
-        np.testing.assert_allclose(
-            edge, np.full(edge.shape, np.nan), equal_nan=True
-        )
-
-
-def test_numpy_equals_dask():
-    small_numpy_based_data_array = xr.DataArray(
-        INPUT_DATA, attrs={'res': (10.0, 10.0)}
-    )
-    small_dask_based_data_array = xr.DataArray(
-        da.from_array(INPUT_DATA, chunks=(2, 2)), attrs={'res': (10.0, 10.0)}
-    )
+    assert_nan_edges_effect(xrspatial_aspect)
 
-    numpy_result = aspect(small_numpy_based_data_array, name='numpy_result')
-    dask_result = aspect(small_dask_based_data_array,
-                         name='dask_result')
-    general_output_checks(small_dask_based_data_array, dask_result)
-    np.testing.assert_allclose(
-        numpy_result.data, dask_result.data.compute(), equal_nan=True)
 
+def test_numpy_equals_dask_qgis_data():
+    # compare using the data run through QGIS
+    numpy_agg = input_data('numpy')
+    dask_agg = input_data('dask+numpy')
+    assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, aspect)
 
-@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_cpu_equals_gpu():
 
-    import cupy
+@pytest.mark.parametrize("size", [(2, 4), (10, 15)])
+@pytest.mark.parametrize(
+    "dtype", [np.int32, np.int64, np.uint32, np.uint64, np.float32, np.float64])
+def test_numpy_equals_dask_random_data(random_data):
+    numpy_agg = create_test_raster(random_data, backend='numpy')
+    dask_agg = create_test_raster(random_data, backend='dask')
+    assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, aspect)
+
+
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+def test_numpy_equals_cupy_qgis_data():
+    # compare using the data run through QGIS
+    numpy_agg = input_data()
+    cupy_agg = input_data('cupy')
+    assert_numpy_equals_cupy(numpy_agg, cupy_agg, aspect)
 
-    small_da = xr.DataArray(INPUT_DATA, attrs={'res': (10.0, 10.0)})
-    small_da_cupy = xr.DataArray(cupy.asarray(INPUT_DATA),
-                                 attrs={'res': (10.0, 10.0)})
 
-    # aspect by xrspatial
-    cpu = aspect(small_da, name='aspect_agg')
-    gpu = aspect(small_da_cupy, name='aspect_agg')
-    general_output_checks(small_da_cupy, gpu)
-    np.testing.assert_allclose(cpu.data, gpu.data.get(), equal_nan=True)
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+@pytest.mark.parametrize("size", [(2, 4), (10, 15)])
+@pytest.mark.parametrize(
+    "dtype", [np.int32, np.int64, np.uint32, np.uint64, np.float32, np.float64])
+def test_numpy_equals_cupy_random_data(random_data):
+    numpy_agg = create_test_raster(random_data, backend='numpy')
+    cupy_agg = create_test_raster(random_data, backend='cupy')
+    assert_numpy_equals_cupy(numpy_agg, cupy_agg, aspect, atol=1e-6, rtol=1e-6)