Repository: makepath/xarray-spatial

Issue:
Improvements to tests
Recently there have been some big improvements to the tests, such as PRs #612 and #613. Here are a few other possible improvements.

  - [x] Tests using 2D arrays should always have `shape[0] != shape[1]`, just in case there are array transpose errors anywhere. e.g. `test_focal.py`s `convolve_2d_data` has shape `(6, 6)`.
  - [x] There are few more `np.isclose(...).all()` that should be replaced with `np.testing.allclose(...)`, e.g. `test_terrain.py` line 39.
  - [ ] There is a `from __future__ import absolute_import` in `test_tiles.py` that I don't think is needed. Such imports were usually related to Python 2/3 changes. But I have not fully researched this.
  - [x] Data for tests should probably use `pytest` fixtures (https://docs.pytest.org/en/6.2.x/fixture.html).

There are some other possible improvements that are really style issues. These affect more code than just the tests, and will probably need some discussion/agreement on policy.

  - [ ] I prefer to group together import of functions from the same file in a single line, rather than a line per import.
  - [x] I would happily increase our max line length from 79 to 100 or so.

Discussion:
thuydotm: Thanks @ianthomas23, I've made the changes as suggested. Can you please run the tests locally to see if it works with rtxpy?

Patch:
diff --git a/benchmarks/asv.conf.json b/benchmarks/asv.conf.json
index 9404957..73d0ad1 100644
--- a/benchmarks/asv.conf.json
+++ b/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.
diff --git a/xrspatial/classify.py b/xrspatial/classify.py
index aa0f47f..36b8df5 100755
--- a/xrspatial/classify.py
+++ b/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
diff --git a/xrspatial/curvature.py b/xrspatial/curvature.py
index 5227141..611048c 100755
--- a/xrspatial/curvature.py
+++ b/xrspatial/curvature.py
@@ -27,7 +27,7 @@ from xrspatial.utils import ArrayTypeFunctionMapping
 
 @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')
diff --git a/xrspatial/proximity.py b/xrspatial/proximity.py
index 15885af..b6f8769 100755
--- a/xrspatial/proximity.py
+++ b/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,
diff --git a/xrspatial/zonal.py b/xrspatial/zonal.py
index dbcbbbb..333749e 100755
--- a/xrspatial/zonal.py
+++ b/xrspatial/zonal.py
@@ -542,7 +542,7 @@ def stats(
         raise ValueError(
             "Got dask-backed DataArray as `values` aggregate. "
             "`stats_funcs` must be a subset of default supported stats "
-            "`[\'mean\', \'max\', \'min\', \'sum\', \'std\', \'var\', \'count\']`"  # noqa
+            "`[\'mean\', \'max\', \'min\', \'sum\', \'std\', \'var\', \'count\']`"
         )
 
     if isinstance(stats_funcs, list):


Test patch:
diff --git a/xrspatial/tests/conftest.py b/xrspatial/tests/conftest.py
new file mode 100644
index 0000000..a56710e
--- /dev/null
+++ b/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
diff --git a/xrspatial/tests/general_checks.py b/xrspatial/tests/general_checks.py
index 377401e..5e97730 100755
--- a/xrspatial/tests/general_checks.py
+++ b/xrspatial/tests/general_checks.py
@@ -2,13 +2,32 @@ import numpy as np
 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
+    )
diff --git a/xrspatial/tests/test_aspect.py b/xrspatial/tests/test_aspect.py
index f3db9cd..413ebdb 100755
--- a/xrspatial/tests/test_aspect.py
+++ b/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)
diff --git a/xrspatial/tests/test_classify.py b/xrspatial/tests/test_classify.py
index b074774..55fc2bb 100755
--- a/xrspatial/tests/test_classify.py
+++ b/xrspatial/tests/test_classify.py
@@ -3,8 +3,6 @@ import pytest
 import xarray as xr
 import numpy as np
 
-import dask.array as da
-
 from xrspatial.utils import doesnt_have_cuda
 from xrspatial import equal_interval
 from xrspatial import natural_breaks
@@ -12,153 +10,126 @@ from xrspatial import quantile
 from xrspatial import reclassify
 
 from xrspatial.tests.general_checks import general_output_checks
+from xrspatial.tests.general_checks import create_test_raster
 
 
-elevation = np.array([
-    [1.,  2.,  3.,  4., np.nan],
-    [5.,  6.,  7.,  8.,  9.],
-    [10., 11., 12., 13., 14.],
-    [15., 16., 17., 18., np.inf],
-])
-
-h, w = elevation.shape
-ys = np.arange(h)
-xs = np.arange(w)
-
-numpy_agg = xr.DataArray(
-    elevation, dims=['y', 'x'], attrs={'res': (10.0, 10.0)}
-)
-numpy_agg['y'] = ys
-numpy_agg['x'] = xs
-
-dask_numpy_agg = xr.DataArray(
-    da.from_array(elevation, chunks=(3, 3)),
-    dims=['y', 'x'],
-    attrs={'res': (10.0, 10.0)}
-)
-dask_numpy_agg['y'] = ys
-dask_numpy_agg['x'] = xs
+def input_data(backend='numpy'):
+    elevation = np.array([
+        [1.,  2.,  3.,  4., np.nan],
+        [5.,  6.,  7.,  8.,  9.],
+        [10., 11., 12., 13., 14.],
+        [15., 16., 17., 18., np.inf],
+    ])
+    raster = create_test_raster(elevation, backend, attrs={'res': (10.0, 10.0)})
+    return raster
 
 
-def test_reclassify_cpu():
+@pytest.fixture
+def result_reclassify():
     bins = [10, 15, np.inf]
     new_values = [1, 2, 3]
-    expected_results = np.asarray([
+    expected_result = np.asarray([
         [1., 1., 1., 1., np.nan],
         [1., 1., 1., 1., 1.],
         [1., 2., 2., 2., 2.],
         [2., 3., 3., 3., 3.]
-    ])
-    # numpy
-    numpy_reclassify = reclassify(
-        numpy_agg, bins=bins, new_values=new_values, name='numpy_reclassify'
-    )
-    general_output_checks(numpy_agg, numpy_reclassify, expected_results)
+    ], dtype=np.float32)
+    return bins, new_values, expected_result
 
-    # dask + numpy
-    dask_reclassify = reclassify(
-        dask_numpy_agg, bins=bins,
-        new_values=new_values, name='dask_reclassify'
-    )
-    general_output_checks(dask_numpy_agg, dask_reclassify, expected_results)
 
+def test_reclassify_numpy(result_reclassify):
+    bins, new_values, expected_result = result_reclassify
+    numpy_agg = input_data()
+    numpy_result = reclassify(numpy_agg, bins=bins, new_values=new_values)
+    general_output_checks(numpy_agg, numpy_result, expected_result, verify_dtype=True)
 
-@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_reclassify_cpu_equals_gpu():
 
-    import cupy
+def test_reclassify_dask_numpy(result_reclassify):
+    bins, new_values, expected_result = result_reclassify
+    dask_agg = input_data(backend='dask')
+    dask_result = reclassify(dask_agg, bins=bins, new_values=new_values)
+    general_output_checks(dask_agg, dask_result, expected_result, verify_dtype=True)
 
-    bins = [10, 15, np.inf]
-    new_values = [1, 2, 3]
 
-    # vanilla numpy version
-    cpu = reclassify(numpy_agg,
-                     name='numpy_result',
-                     bins=bins,
-                     new_values=new_values)
-
-    # cupy
-    cupy_agg = xr.DataArray(cupy.asarray(elevation),
-                            attrs={'res': (10.0, 10.0)})
-    gpu = reclassify(cupy_agg,
-                     name='cupy_result',
-                     bins=bins,
-                     new_values=new_values)
-    general_output_checks(cupy_agg, gpu)
-    np.testing.assert_allclose(cpu.data, gpu.data.get(), equal_nan=True)
-
-    # dask + cupy
-    dask_cupy_agg = xr.DataArray(
-        cupy.asarray(elevation), attrs={'res': (10.0, 10.0)})
-    dask_cupy_agg.data = da.from_array(dask_cupy_agg.data, chunks=(3, 3))
-    dask_gpu = reclassify(
-        dask_cupy_agg, name='dask_cupy_result',
-        bins=bins, new_values=new_values
-    )
-    general_output_checks(dask_cupy_agg, dask_gpu)
-    dask_gpu.data = dask_gpu.data.compute()
-    np.testing.assert_allclose(cpu.data, dask_gpu.data.get(), equal_nan=True)
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+def test_reclassify_cupy(result_reclassify):
+    bins, new_values, expected_result = result_reclassify
+    cupy_agg = input_data(backend='cupy')
+    cupy_result = reclassify(cupy_agg, bins=bins, new_values=new_values)
+    general_output_checks(cupy_agg, cupy_result, expected_result, verify_dtype=True)
+
+
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+def test_reclassify_dask_cupy(result_reclassify):
+    bins, new_values, expected_result = result_reclassify
+    dask_cupy_agg = input_data(backend='dask+cupy')
+    dask_cupy_result = reclassify(dask_cupy_agg, bins=bins, new_values=new_values)
+    general_output_checks(dask_cupy_agg, dask_cupy_result, expected_result, verify_dtype=True)
 
 
-def test_quantile_cpu():
+@pytest.fixture
+def result_quantile():
     k = 5
-    expected_results = np.asarray([
+    expected_result = np.asarray([
         [0., 0., 0., 0., np.nan],
         [1., 1., 1., 2., 2.],
         [2., 2., 3., 3., 3.],
         [4., 4., 4., 4., np.nan]
-    ])
-    # numpy
+    ], dtype=np.float32)
+    return k, expected_result
+
+
+def test_quantile_numpy(result_quantile):
+    k, expected_result = result_quantile
+    numpy_agg = input_data()
     numpy_quantile = quantile(numpy_agg, k=k)
-    general_output_checks(numpy_agg, numpy_quantile, expected_results)
+    general_output_checks(numpy_agg, numpy_quantile, expected_result, verify_dtype=True)
 
-    # dask + numpy
-    dask_quantile = quantile(dask_numpy_agg, k=k)
-    general_output_checks(dask_numpy_agg, dask_quantile)
-    dask_quantile = dask_quantile.compute()
-    unique_elements = np.unique(
-        dask_quantile.data[np.isfinite(dask_quantile.data)]
-    )
-    assert len(unique_elements) == k
 
+def test_quantile_dask_numpy(result_quantile):
     #     Note that dask's percentile algorithm is
     #     approximate, while numpy's is exact.
     #     This may cause some differences between
     #     results of vanilla numpy and
     #     dask version of the input agg.
     #     https://github.com/dask/dask/issues/3099
-    #     This assertion may fail
-    # dask_quantile = dask_quantile.compute()
-    # assert np.isclose(numpy_quantile, dask_quantile, equal_nan=True).all()
-
 
-@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_quantile_cpu_equals_gpu():
+    dask_numpy_agg = input_data('dask+numpy')
+    k, expected_result = result_quantile
+    dask_quantile = quantile(dask_numpy_agg, k=k)
+    general_output_checks(dask_numpy_agg, dask_quantile)
+    dask_quantile = dask_quantile.compute()
+    unique_elements = np.unique(
+        dask_quantile.data[np.isfinite(dask_quantile.data)]
+    )
+    assert len(unique_elements) == k
 
-    import cupy
 
-    k = 5
-    # vanilla numpy version
-    cpu = quantile(numpy_agg, k=k, name='numpy_result')
-    # cupy
-    cupy_agg = xr.DataArray(
-        cupy.asarray(elevation), attrs={'res': (10.0, 10.0)})
-    gpu = quantile(cupy_agg, k=k, name='cupy_result')
-    general_output_checks(cupy_agg, gpu)
-    np.testing.assert_allclose(cpu.data, gpu.data.get(), equal_nan=True)
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+def test_quantile_cupy(result_quantile):
+    k, expected_result = result_quantile
+    cupy_agg = input_data('cupy')
+    cupy_result = quantile(cupy_agg, k=k)
+    general_output_checks(cupy_agg, cupy_result, expected_result, verify_dtype=True)
 
 
-def test_natural_breaks_cpu():
+@pytest.fixture
+def result_natural_breaks():
     k = 5
-    expected_results = np.asarray([
+    expected_result = np.asarray([
         [0., 0., 0., 1., np.nan],
         [1., 1., 2., 2., 2.],
         [2., 3., 3., 3., 3.],
         [4., 4., 4., 4., np.nan]
-    ])
-    # vanilla numpy
+    ], dtype=np.float32)
+    return k, expected_result
+
+
+def test_natural_breaks_numpy(result_natural_breaks):
+    numpy_agg = input_data()
+    k, expected_result = result_natural_breaks
     numpy_natural_breaks = natural_breaks(numpy_agg, k=k)
-    general_output_checks(numpy_agg, numpy_natural_breaks, expected_results)
+    general_output_checks(numpy_agg, numpy_natural_breaks, expected_result, verify_dtype=True)
 
 
 def test_natural_breaks_cpu_deterministic():
@@ -186,52 +157,42 @@ def test_natural_breaks_cpu_deterministic():
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_natural_breaks_cpu_equals_gpu():
-
-    import cupy
+def test_natural_breaks_cupy(result_natural_breaks):
+    cupy_agg = input_data('cupy')
+    k, expected_result = result_natural_breaks
+    cupy_natural_breaks = natural_breaks(cupy_agg, k=k)
+    general_output_checks(cupy_agg, cupy_natural_breaks, expected_result, verify_dtype=True)
 
-    k = 5
-    # vanilla numpy version
-    cpu = natural_breaks(numpy_agg, k=k, name='numpy_result')
-    # cupy
-    cupy_agg = xr.DataArray(cupy.asarray(elevation),
-                            attrs={'res': (10.0, 10.0)})
-    gpu = natural_breaks(cupy_agg, k=k, name='cupy_result')
-    general_output_checks(cupy_agg, gpu)
-    np.testing.assert_allclose(cpu.data, gpu.data.get(), equal_nan=True)
 
-
-def test_equal_interval_cpu():
+@pytest.fixture
+def result_equal_interval():
     k = 3
-    expected_results = np.asarray([
+    expected_result = np.asarray([
         [0., 0., 0., 0., np.nan],
         [0., 0., 1., 1., 1.],
         [1., 1., 1., 2., 2.],
         [2., 2., 2., 2., np.nan]
-    ])
+    ], dtype=np.float32)
+    return k, expected_result
 
-    # numpy
-    numpy_ei = equal_interval(numpy_agg, k=k)
-    general_output_checks(numpy_agg, numpy_ei, expected_results)
 
-    # dask + numpy
-    dask_ei = equal_interval(dask_numpy_agg, k=k, name='dask_reclassify')
-    general_output_checks(dask_numpy_agg, dask_ei)
-    dask_ei = dask_ei.compute()
-    np.testing.assert_allclose(numpy_ei, dask_ei, equal_nan=True)
+def test_equal_interval_numpy(result_equal_interval):
+    k, expected_result = result_equal_interval
+    numpy_agg = input_data('numpy')
+    numpy_result = equal_interval(numpy_agg, k=k)
+    general_output_checks(numpy_agg, numpy_result, expected_result, verify_dtype=True)
 
 
-@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_equal_interval_cpu_equals_gpu():
+def test_equal_interval_dask_numpy(result_equal_interval):
+    k, expected_result = result_equal_interval
+    dask_agg = input_data('dask+numpy')
+    dask_numpy_result = equal_interval(dask_agg, k=k)
+    general_output_checks(dask_agg, dask_numpy_result, expected_result, verify_dtype=True)
 
-    import cupy
 
-    k = 3
-    # numpy
-    cpu = equal_interval(numpy_agg, k=k)
-    # cupy
-    cupy_agg = xr.DataArray(
-        cupy.asarray(elevation), attrs={'res': (10.0, 10.0)})
-    gpu = equal_interval(cupy_agg, k=k)
-    general_output_checks(cupy_agg, gpu)
-    np.testing.assert_allclose(cpu.data, gpu.data.get(), equal_nan=True)
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+def test_equal_interval_cupy(result_equal_interval):
+    k, expected_result = result_equal_interval
+    cupy_agg = input_data(backend='cupy')
+    cupy_result = equal_interval(cupy_agg, k=k)
+    general_output_checks(cupy_agg, cupy_result, expected_result, verify_dtype=True)
diff --git a/xrspatial/tests/test_curvature.py b/xrspatial/tests/test_curvature.py
index 1e373ae..c019d22 100755
--- a/xrspatial/tests/test_curvature.py
+++ b/xrspatial/tests/test_curvature.py
@@ -1,115 +1,104 @@
 import pytest
 import numpy as np
-import xarray as xr
-
-import dask.array as da
 
 from xrspatial import curvature
 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 general_output_checks
 
 
-elevation = 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
-)
-
-
-def test_curvature_on_flat_surface():
-    # flat surface
-    test_arr1 = np.array([[0, 0, 0, 0, 0],
-                          [0, 0, 0, 0, 0],
-                          [0, 0, 0, 0, 0],
-                          [0, 0, 0, 0, 0],
-                          [0, 0, 0, 0, 0]])
-    expected_results = np.array([
-        [np.nan, np.nan, np.nan, np.nan, np.nan],
-        [np.nan,      0,      0,      0, np.nan],
-        [np.nan,      0,      0,      0, np.nan],
-        [np.nan,      0,      0,      0, np.nan],
-        [np.nan, np.nan, np.nan, np.nan, np.nan]
-    ])
-    test_raster1 = xr.DataArray(test_arr1, attrs={'res': (1, 1)})
-    curv = curvature(test_raster1)
-    general_output_checks(test_raster1, curv, expected_results)
-
-
-def test_curvature_on_convex_surface():
-    # convex
-    test_arr2 = np.array([[0, 0, 0, 0, 0],
-                          [0, 0, 0, 0, 0],
-                          [0, 0, -1, 0, 0],
-                          [0, 0, 0, 0, 0],
-                          [0, 0, 0, 0, 0]])
-    expected_results = np.asarray([
-         [np.nan, np.nan, np.nan,  np.nan, np.nan],
-         [np.nan,     0.,   100.,      0., np.nan],
-         [np.nan,   100.,  -400.,    100., np.nan],
-         [np.nan,     0.,   100.,      0., np.nan],
-         [np.nan,  np.nan, np.nan, np.nan, np.nan]
-    ])
-    test_raster2 = xr.DataArray(test_arr2, attrs={'res': (1, 1)})
-    curv = curvature(test_raster2)
-    general_output_checks(test_raster2, curv, expected_results)
-
-
-def test_curvature_on_concave_surface():
-    # concave
-    test_arr3 = np.array([[0, 0, 0, 0, 0],
-                          [0, 0, 0, 0, 0],
-                          [0, 0, 1, 0, 0],
-                          [0, 0, 0, 0, 0],
-                          [0, 0, 0, 0, 0]])
-    expected_results = np.asarray([
-         [np.nan, np.nan, np.nan,  np.nan, np.nan],
-         [np.nan,     0.,  -100.,      0., np.nan],
-         [np.nan,  -100.,   400.,   -100., np.nan],
-         [np.nan,     0.,  -100.,      0., np.nan],
-         [np.nan,  np.nan, np.nan, np.nan, np.nan]
-    ])
-    test_raster3 = xr.DataArray(test_arr3, attrs={'res': (1, 1)})
-    curv = curvature(test_raster3)
-    general_output_checks(test_raster3, curv, expected_results)
+@pytest.fixture
+def flat_surface(size, dtype):
+    flat = np.zeros(size, dtype=dtype)
+    expected_result = np.zeros(size, dtype=np.float32)
+    # nan edges effect
+    expected_result[0, :] = np.nan
+    expected_result[-1, :] = np.nan
+    expected_result[:, 0] = np.nan
+    expected_result[:, -1] = np.nan
+    return flat, expected_result
+
+
+@pytest.fixture
+def convex_surface():
+    convex_data = np.array([
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 0, -1, 0, 0],
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0]])
+    expected_result = np.asarray([
+         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
+         [np.nan, 0,      0.,     100.,     0.,   np.nan],
+         [np.nan, 0,      100.,  -400.,   100.,   np.nan],
+         [np.nan, 0,      0.,     100.,     0.,   np.nan],
+         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]
+    ], dtype=np.float32)
+    return convex_data, expected_result
+
+
+@pytest.fixture
+def concave_surface():
+    concave_data = np.array([
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 1, 0, 0],
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0]])
+    expected_result = np.asarray([
+         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
+         [np.nan, 0,      0.,    -100.,     0.,   np.nan],
+         [np.nan, 0,     -100.,   400.,  -100.,   np.nan],
+         [np.nan, 0,      0.,    -100.,     0.,   np.nan],
+         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]
+    ], dtype=np.float32)
+    return concave_data, expected_result
+
+
+@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_curvature_on_flat_surface(flat_surface):
+    flat_data, expected_result = flat_surface
+    numpy_agg = create_test_raster(flat_data, attrs={'res': (1, 1)})
+    numpy_result = curvature(numpy_agg)
+    general_output_checks(numpy_agg, numpy_result, expected_result, verify_dtype=True)
+
+
+def test_curvature_on_convex_surface(convex_surface):
+    convex_data, expected_result = convex_surface
+    numpy_agg = create_test_raster(convex_data, attrs={'res': (1, 1)})
+    numpy_result = curvature(numpy_agg)
+    general_output_checks(numpy_agg, numpy_result, expected_result, verify_dtype=True)
+
+
+def test_curvature_on_concave_surface(concave_surface):
+    concave_data, expected_result = concave_surface
+    numpy_agg = create_test_raster(concave_data, attrs={'res': (1, 1)})
+    numpy_result = curvature(numpy_agg)
+    general_output_checks(numpy_agg, numpy_result, expected_result, verify_dtype=True)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_curvature_gpu_equals_cpu():
-
-    import cupy
-
-    agg_numpy = xr.DataArray(elevation, attrs={'res': (10.0, 10.0)})
-    cpu = curvature(agg_numpy, name='numpy_result')
-
-    agg_cupy = xr.DataArray(
-        cupy.asarray(elevation), attrs={'res': (10.0, 10.0)}
-    )
-    gpu = curvature(agg_cupy, name='cupy_result')
-
-    general_output_checks(agg_cupy, gpu)
-    np.testing.assert_allclose(cpu.data, gpu.data.get(), equal_nan=True)
-
+@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, curvature)
     # NOTE: Dask + GPU code paths don't currently work because of
     # dask casting cupy arrays to numpy arrays during
     # https://github.com/dask/dask/issues/4842
 
 
-def test_curvature_numpy_equals_dask():
-    agg_numpy = xr.DataArray(elevation, attrs={'res': (10.0, 10.0)})
-    numpy_curvature = curvature(agg_numpy, name='numpy_curvature')
-
-    agg_dask = xr.DataArray(
-        da.from_array(elevation, chunks=(3, 3)), attrs={'res': (10.0, 10.0)}
-    )
-    dask_curvature = curvature(agg_dask, name='dask_curvature')
-    general_output_checks(agg_dask, dask_curvature)
-
-    # both produce same results
-    np.testing.assert_allclose(
-        numpy_curvature.data, dask_curvature.data.compute(), equal_nan=True)
+@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, curvature)
diff --git a/xrspatial/tests/test_focal.py b/xrspatial/tests/test_focal.py
index a2a4e0c..ef9b2af 100755
--- a/xrspatial/tests/test_focal.py
+++ b/xrspatial/tests/test_focal.py
@@ -10,9 +10,10 @@ from xrspatial.utils import ngjit
 from xrspatial import mean
 from xrspatial.focal import hotspots, apply, focal_stats
 from xrspatial.convolution import (
-    convolve_2d, calc_cellsize, circle_kernel, annulus_kernel
+    convolve_2d, circle_kernel, annulus_kernel
 )
 
+from xrspatial.tests.general_checks import create_test_raster
 from xrspatial.tests.general_checks import general_output_checks
 
 
@@ -88,65 +89,34 @@ def test_mean_transfer_function_gpu_equals_cpu():
         assert e_info
 
 
-convolve_2d_data = np.array([[0., 1., 1., 1., 1., 1.],
-                             [1., 0., 1., 1., 1., 1.],
-                             [1., 1., 0., 1., 1., 1.],
-                             [1., 1., 1., np.nan, 1., 1.],
-                             [1., 1., 1., 1., 0., 1.],
-                             [1., 1., 1., 1., 1., 0.]])
-
-
-def test_kernel():
-    data = convolve_2d_data
-    m, n = data.shape
-    agg = xr.DataArray(data, dims=['y', 'x'])
-    agg['x'] = np.linspace(0, n, n)
-    agg['y'] = np.linspace(0, m, m)
-
-    cellsize_x, cellsize_y = calc_cellsize(agg)
-
-    kernel1 = circle_kernel(cellsize_x, cellsize_y, 2)
-    expected_kernel1 = np.array([[0, 1, 0],
-                                 [1, 1, 1],
-                                 [0, 1, 0]])
-    assert isinstance(kernel1, np.ndarray)
-    np.testing.assert_allclose(kernel1, expected_kernel1, equal_nan=True)
-
-    kernel2 = annulus_kernel(cellsize_x, cellsize_y, 2, 0.5)
-    expected_kernel2 = np.array([[0, 1, 0],
-                                 [1, 0, 1],
-                                 [0, 1, 0]])
-    assert isinstance(kernel2, np.ndarray)
-    np.testing.assert_allclose(kernel2, expected_kernel2, equal_nan=True)
-
-
-def test_convolution():
-    data = convolve_2d_data
-    dask_data = da.from_array(data, chunks=(3, 3))
-
-    kernel1 = np.ones((1, 1))
-    numpy_output_1 = convolve_2d(data, kernel1)
-    expected_output_1 = np.array([[0., 1., 1., 1., 1., 1.],
-                                  [1., 0., 1., 1., 1., 1.],
-                                  [1., 1., 0., 1., 1., 1.],
-                                  [1., 1., 1., np.nan, 1., 1.],
-                                  [1., 1., 1., 1., 0., 1.],
-                                  [1., 1., 1., 1., 1., 0.]])
-    assert isinstance(numpy_output_1, np.ndarray)
-    np.testing.assert_allclose(
-        numpy_output_1, expected_output_1, equal_nan=True)
+@pytest.fixture
+def convolve_2d_data():
+    data = np.array([
+        [0., 1., 1., 1., 1., 1.],
+        [1., 0., 1., 1., 1., 1.],
+        [1., 1., 0., 1., 1., 1.],
+        [1., 1., 1., np.nan, 1., 1.],
+        [1., 1., 1., 1., 0., 1.],
+        [1., 1., 1., 1., 1., 0.]
+    ])
+    return data
 
-    dask_output_1 = convolve_2d(dask_data, kernel1)
-    assert isinstance(dask_output_1, da.Array)
-    np.testing.assert_allclose(
-        dask_output_1.compute(), expected_output_1, equal_nan=True
-    )
 
-    kernel2 = np.array([[0, 1, 0],
-                        [1, 1, 1],
-                        [0, 1, 0]])
-    numpy_output_2 = convolve_2d(data, kernel2)
-    expected_output_2 = np.array([
+@pytest.fixture
+def kernel_circle_1_1_1():
+    result = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
+    return result
+
+
+@pytest.fixture
+def kernel_annulus_2_2_2_1():
+    result = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]])
+    return result
+
+
+@pytest.fixture
+def convolution_kernel_circle_1_1_1():
+    expected_result = np.array([
         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
         [np.nan, 4., 3., 5., 5., np.nan],
         [np.nan, 3., np.nan, np.nan, np.nan, np.nan],
@@ -154,24 +124,12 @@ def test_convolution():
         [np.nan, 5., np.nan, np.nan, np.nan, np.nan],
         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]
     ])
-    # kernel2 is of 3x3, thus the border edge is 1 cell long.
-    # currently, ignoring border edge (i.e values in edges are all nans)
-    assert isinstance(numpy_output_2, np.ndarray)
-    np.testing.assert_allclose(
-        numpy_output_2, expected_output_2, equal_nan=True
-    )
+    return expected_result
 
-    dask_output_2 = convolve_2d(dask_data, kernel2)
-    assert isinstance(dask_output_2, da.Array)
-    np.testing.assert_allclose(
-        dask_output_2.compute(), expected_output_2, equal_nan=True
-    )
 
-    kernel3 = np.array([[0, 1, 0],
-                        [1, 0, 1],
-                        [0, 1, 0]])
-    numpy_output_3 = convolve_2d(data, kernel3)
-    expected_output_3 = np.array([
+@pytest.fixture
+def convolution_kernel_annulus_2_2_1():
+    expected_result = np.array([
         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
         [np.nan, 4., 2., 4., 4., np.nan],
         [np.nan, 2., np.nan, np.nan, np.nan, np.nan],
@@ -179,125 +137,169 @@ def test_convolution():
         [np.nan, 4., np.nan, np.nan, np.nan, np.nan],
         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]
     ])
-    # kernel3 is of 3x3, thus the border edge is 1 cell long.
-    # currently, values in edges are all nans
-    assert isinstance(numpy_output_3, np.ndarray)
+    return expected_result
+
+
+def test_kernel(kernel_circle_1_1_1, kernel_annulus_2_2_2_1):
+    kernel_circle = circle_kernel(1, 1, 1)
+    assert isinstance(kernel_circle, np.ndarray)
+    np.testing.assert_allclose(kernel_circle, kernel_circle_1_1_1, equal_nan=True)
+
+    kernel_annulus = annulus_kernel(2, 2, 2, 1)
+    assert isinstance(kernel_annulus, np.ndarray)
+    np.testing.assert_allclose(kernel_annulus, kernel_annulus_2_2_2_1, equal_nan=True)
+
+
+def test_convolution_numpy(
+    convolve_2d_data,
+    kernel_circle_1_1_1,
+    convolution_kernel_circle_1_1_1,
+    kernel_annulus_2_2_2_1,
+    convolution_kernel_annulus_2_2_1
+):
+    kernel_custom = np.ones((1, 1))
+    result_kernel_custom = convolve_2d(convolve_2d_data, kernel_custom)
+    assert isinstance(result_kernel_custom, np.ndarray)
+    # kernel is [[1]], thus the result equals input data
+    np.testing.assert_allclose(result_kernel_custom, convolve_2d_data, equal_nan=True)
+
+    result_kernel_circle = convolve_2d(convolve_2d_data, kernel_circle_1_1_1)
+    assert isinstance(result_kernel_circle, np.ndarray)
     np.testing.assert_allclose(
-        numpy_output_3, expected_output_3, equal_nan=True)
+        result_kernel_circle, convolution_kernel_circle_1_1_1, equal_nan=True
+    )
 
-    dask_output_3 = convolve_2d(dask_data, kernel3)
-    assert isinstance(dask_output_3, da.Array)
+    result_kernel_annulus = convolve_2d(convolve_2d_data, kernel_annulus_2_2_2_1)
+    assert isinstance(result_kernel_annulus, np.ndarray)
     np.testing.assert_allclose(
-        dask_output_3.compute(), expected_output_3, equal_nan=True
+        result_kernel_annulus, convolution_kernel_annulus_2_2_1, equal_nan=True
     )
 
 
-@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_2d_convolution_gpu_equals_cpu():
+def test_convolution_dask_numpy(
+    convolve_2d_data,
+    kernel_circle_1_1_1,
+    convolution_kernel_circle_1_1_1,
+    kernel_annulus_2_2_2_1,
+    convolution_kernel_annulus_2_2_1
+):
+    dask_data = da.from_array(convolve_2d_data, chunks=(3, 3))
+    kernel_custom = np.ones((1, 1))
+    result_kernel_custom = convolve_2d(dask_data, kernel_custom)
+    assert isinstance(result_kernel_custom, da.Array)
+    # kernel is [[1]], thus the result equals input data
+    np.testing.assert_allclose(result_kernel_custom.compute(), convolve_2d_data, equal_nan=True)
+
+    result_kernel_circle = convolve_2d(dask_data, kernel_circle_1_1_1)
+    assert isinstance(result_kernel_circle, da.Array)
+    np.testing.assert_allclose(
+        result_kernel_circle.compute(), convolution_kernel_circle_1_1_1, equal_nan=True
+    )
 
+    result_kernel_annulus = convolve_2d(dask_data, kernel_annulus_2_2_2_1)
+    assert isinstance(result_kernel_annulus, da.Array)
+    np.testing.assert_allclose(
+        result_kernel_annulus.compute(), convolution_kernel_annulus_2_2_1, equal_nan=True
+    )
+
+
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+def test_2d_convolution_gpu(
+    convolve_2d_data,
+    kernel_circle_1_1_1,
+    convolution_kernel_circle_1_1_1,
+    kernel_annulus_2_2_2_1,
+    convolution_kernel_annulus_2_2_1
+):
     import cupy
+    cupy_data = cupy.asarray(convolve_2d_data)
 
-    data = convolve_2d_data
-    numpy_agg = xr.DataArray(data)
-    cupy_agg = xr.DataArray(cupy.asarray(data))
+    kernel_custom = np.ones((1, 1))
+    result_kernel_custom = convolve_2d(cupy_data, kernel_custom)
+    assert isinstance(result_kernel_custom, cupy.ndarray)
+    # kernel is [[1]], thus the result equals input data
+    np.testing.assert_allclose(result_kernel_custom.get(), convolve_2d_data, equal_nan=True)
 
-    kernel1 = np.ones((1, 1))
-    output_numpy1 = convolve_2d(numpy_agg.data, kernel1)
-    output_cupy1 = convolve_2d(cupy_agg.data, kernel1)
-    assert isinstance(output_cupy1, cupy.ndarray)
+    result_kernel_circle = convolve_2d(cupy_data, kernel_circle_1_1_1)
+    assert isinstance(result_kernel_circle, cupy.ndarray)
     np.testing.assert_allclose(
-        output_numpy1, output_cupy1.get(), equal_nan=True)
-
-    kernel2 = np.array([[0, 1, 0],
-                        [1, 1, 1],
-                        [0, 1, 0]])
-    output_numpy2 = convolve_2d(numpy_agg.data, kernel2)
-    output_cupy2 = convolve_2d(cupy_agg.data, kernel2)
-    assert isinstance(output_cupy2, cupy.ndarray)
-    np.testing.assert_allclose(
-        output_numpy2, output_cupy2.get(), equal_nan=True)
-
-    kernel3 = np.array([[0, 1, 0],
-                        [1, 0, 1],
-                        [0, 1, 0]])
-    output_numpy3 = convolve_2d(numpy_agg.data, kernel3)
-    output_cupy3 = convolve_2d(cupy_agg.data, kernel3)
-    assert isinstance(output_cupy3, cupy.ndarray)
+        result_kernel_circle.get(), convolution_kernel_circle_1_1_1, equal_nan=True
+    )
+
+    result_kernel_annulus = convolve_2d(cupy_data, kernel_annulus_2_2_2_1)
+    assert isinstance(result_kernel_annulus, cupy.ndarray)
     np.testing.assert_allclose(
-        output_numpy3, output_cupy3.get(), equal_nan=True)
+        result_kernel_annulus.get(), convolution_kernel_annulus_2_2_1, equal_nan=True
+    )
 
     # dask + cupy case not implemented
     dask_cupy_agg = xr.DataArray(
-        da.from_array(cupy.asarray(data), chunks=(3, 3))
+        da.from_array(cupy.asarray(convolve_2d_data), chunks=(3, 3))
     )
     with pytest.raises(NotImplementedError) as e_info:
-        convolve_2d(dask_cupy_agg.data, kernel3)
+        convolve_2d(dask_cupy_agg.data, kernel_custom)
         assert e_info
 
 
-data_apply = np.array([[0, 1, 2, 3, 4, 5],
-                       [6, 7, 8, 9, 10, 11],
-                       [12, 13, 14, 15, 16, 17],
-                       [18, 19, 20, 21, 22, 23]])
-kernel_apply = np.array([
-    [0, 1, 0],
-    [1, 0, 1],
-    [0, 1, 0]
-])
+@pytest.fixture
+def data_apply():
+    data = np.array([[0, 1, 2, 3, 4, 5],
+                     [6, 7, 8, 9, 10, 11],
+                     [12, 13, 14, 15, 16, 17],
+                     [18, 19, 20, 21, 22, 23]])
+    kernel = np.array([
+        [0, 1, 0],
+        [1, 0, 1],
+        [0, 1, 0]
+    ])
+    expected_result = np.zeros_like(data)
+    return data, kernel, expected_result
 
 
-def test_apply_cpu():
-    @ngjit
-    def func_zero_cpu(x):
-        return 0
-    expected_results = np.zeros_like(data_apply)
+def func_zero(x):
+    return 0
 
-    # numpy case
-    numpy_agg = xr.DataArray(data_apply)
-    numpy_apply = apply(numpy_agg, kernel_apply, func_zero_cpu)
-    general_output_checks(numpy_agg, numpy_apply, expected_results)
 
-    # dask + numpy case
-    dask_numpy_agg = xr.DataArray(da.from_array(data_apply, chunks=(3, 3)))
-    dask_numpy_apply = apply(dask_numpy_agg, kernel_apply, func_zero_cpu)
-    general_output_checks(dask_numpy_agg, dask_numpy_apply, expected_results)
+@ngjit
+def func_zero_cpu(x):
+    return 0
 
 
-@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_apply_gpu():
-    def func_zero(x):
-        return 0
+def test_apply_numpy(data_apply):
+    data, kernel, expected_result = data_apply
+    numpy_agg = create_test_raster(data)
+    numpy_apply = apply(numpy_agg, kernel, func_zero_cpu)
+    general_output_checks(numpy_agg, numpy_apply, expected_result)
 
-    @ngjit
-    def func_zero_cpu(x):
-        return 0
 
-    expected_results = np.zeros_like(data_apply)
+def test_apply_dask_numpy(data_apply):
+    data, kernel, expected_result = data_apply
+    dask_numpy_agg = create_test_raster(data, backend='dask')
+    dask_numpy_apply = apply(dask_numpy_agg, kernel, func_zero_cpu)
+    general_output_checks(dask_numpy_agg, dask_numpy_apply, expected_result)
 
+
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+def test_apply_gpu(data_apply):
+    data, kernel, expected_result = data_apply
     # cupy case
-    import cupy
-    cupy_agg = xr.DataArray(cupy.asarray(data_apply))
-    cupy_apply = apply(cupy_agg, kernel_apply, func_zero)
-    general_output_checks(cupy_agg, cupy_apply, expected_results)
+    cupy_agg = create_test_raster(data, backend='cupy')
+    cupy_apply = apply(cupy_agg, kernel, func_zero)
+    general_output_checks(cupy_agg, cupy_apply, expected_result)
 
     # dask + cupy case not implemented
-    dask_cupy_agg = xr.DataArray(
-        da.from_array(cupy.asarray(data_apply), chunks=(3, 3))
-    )
+    dask_cupy_agg = create_test_raster(data, backend='dask+cupy')
     with pytest.raises(NotImplementedError) as e_info:
-        apply(dask_cupy_agg, kernel_apply, func_zero)
+        apply(dask_cupy_agg, kernel, func_zero)
         assert e_info
 
 
-def test_focal_stats_cpu():
+@pytest.fixture
+def data_focal_stats():
     data = np.arange(16).reshape(4, 4)
-    numpy_agg = xr.DataArray(data)
-    dask_numpy_agg = xr.DataArray(da.from_array(data, chunks=(3, 3)))
-
     cellsize = (1, 1)
     kernel = circle_kernel(*cellsize, 1.5)
-
-    expected_results = np.asarray([
+    expected_result = np.asarray([
         # mean
         [[1.66666667, 2., 3., 4.],
          [4.25, 5., 6., 6.75],
@@ -334,23 +336,30 @@ def test_focal_stats_cpu():
          [33., 45., 50., 43.],
          [33., 48., 52., 40.]]
     ])
+    return data, kernel, expected_result
 
+
+def test_focal_stats_numpy(data_focal_stats):
+    data, kernel, expected_result = data_focal_stats
+    numpy_agg = create_test_raster(data)
     numpy_focalstats = focal_stats(numpy_agg, kernel)
     general_output_checks(
-        numpy_agg, numpy_focalstats,
-        verify_attrs=False, expected_results=expected_results
+        numpy_agg, numpy_focalstats, verify_attrs=False, expected_results=expected_result
     )
     assert numpy_focalstats.ndim == 3
 
+
+def test_focal_stats_dask_numpy(data_focal_stats):
+    data, kernel, expected_result = data_focal_stats
+    dask_numpy_agg = create_test_raster(data, backend='dask')
     dask_numpy_focalstats = focal_stats(dask_numpy_agg, kernel)
     general_output_checks(
-        dask_numpy_agg, dask_numpy_focalstats,
-        verify_attrs=False, expected_results=expected_results
+        dask_numpy_agg, dask_numpy_focalstats, verify_attrs=False, expected_results=expected_result
     )
 
 
-def test_hotspots():
-    n, m = 10, 10
+@pytest.fixture
+def data_hotspots():
     data = np.asarray([
         [np.nan, 0., 0., 0., 0., 0., 0., 0., 0., 0.],
         [0., 10000., 10000., 10000., 0., 0., 0., 0., 0., 0.],
@@ -363,74 +372,46 @@ def test_hotspots():
         [0., 0., 0., 0., 0., 0., 0., -10000., -10000., -10000.],
         [0., 0., 0., 0., 0., 0., 0., -10000., -10000., -10000.]
     ])
-    numpy_agg = xr.DataArray(data, dims=['y', 'x'])
-    numpy_agg['x'] = np.linspace(0, n, n)
-    numpy_agg['y'] = np.linspace(0, m, m)
-
-    cellsize_x, cellsize_y = calc_cellsize(numpy_agg)
-    kernel = circle_kernel(cellsize_x, cellsize_y, 2.0)
-
-    expected_results = np.array([
-        [0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
-        [0,  0, 90,  0,  0,  0,  0,  0,  0,  0],
-        [0, 90, 95, 90,  0,  0,  0,  0,  0,  0],
-        [0,  0, 90,  0,  0,  0,  0,  0,  0,  0],
-        [0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
-        [0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
-        [0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
-        [0,  0,  0,  0,  0,  0,  0,  0, -90,  0],
-        [0,  0,  0,  0,  0,  0,  0, -90, -95,  0],
-        [0,  0,  0,  0,  0,  0,  0,  0,  0,  0]
+    kernel = np.array([[0., 1., 0.], [1., 1., 1.], [0., 1., 0.]])
+    expected_result = np.array([
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 90, 0, 0, 0, 0, 0, 0, 0],
+        [0, 90, 95, 90, 0, 0, 0, 0, 0, 0],
+        [0, 0, 90, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, -90, 0],
+        [0, 0, 0, 0, 0, 0, 0, -90, -95, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
     ], dtype=np.int8)
 
-    # numpy case
-    numpy_hotspots = hotspots(numpy_agg, kernel)
-    general_output_checks(numpy_agg, numpy_hotspots, expected_results)
+    return data, kernel, expected_result
 
-    # dask + numpy
-    dask_numpy_agg = xr.DataArray(da.from_array(data, chunks=(3, 3)))
-    dask_numpy_hotspots = hotspots(dask_numpy_agg, kernel)
-    general_output_checks(
-        dask_numpy_agg, dask_numpy_hotspots, expected_results)
 
+def test_hotspots_numpy(data_hotspots):
+    data, kernel, expected_result = data_hotspots
+    numpy_agg = create_test_raster(data)
+    numpy_hotspots = hotspots(numpy_agg, kernel)
+    general_output_checks(numpy_agg, numpy_hotspots, expected_result)
 
-@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_hotspot_gpu_equals_cpu():
-    n, m = 10, 10
-    data = np.asarray([
-        [np.nan, 0., 0., 0., 0., 0., 0., 0., 0., 0.],
-        [0., 10000., 10000., 10000., 0., 0., 0., 0., 0., 0.],
-        [0., 10000., 10000., 10000., 0., 0., 0., 0., 0., 0.],
-        [0., 10000., 10000., 10000., 0., 0., 0., 0., 0., 0.],
-        [0., 0., 0., 0., np.nan, 0., 0., 0., 0., 0.],
-        [0., 0., 0., 0., 0., np.nan, 0., 0., 0., 0.],
-        [0., 0., 0., 0., 0., 0., np.nan, 0., 0., 0.],
-        [0., 0., 0., 0., 0., 0., 0., -10000., -10000., -10000.],
-        [0., 0., 0., 0., 0., 0., 0., -10000., -10000., -10000.],
-        [0., 0., 0., 0., 0., 0., 0., -10000., -10000., -10000.]
-    ])
-    numpy_agg = xr.DataArray(data, dims=['y', 'x'])
-    numpy_agg['x'] = np.linspace(0, n, n)
-    numpy_agg['y'] = np.linspace(0, m, m)
 
-    cellsize_x, cellsize_y = calc_cellsize(numpy_agg)
-    kernel = circle_kernel(cellsize_x, cellsize_y, 2.0)
+def test_hotspots_dask_numpy(data_hotspots):
+    data, kernel, expected_result = data_hotspots
+    dask_numpy_agg = create_test_raster(data, backend='dask')
+    dask_numpy_hotspots = hotspots(dask_numpy_agg, kernel)
+    general_output_checks(dask_numpy_agg, dask_numpy_hotspots, expected_result)
 
-    # numpy case
-    numpy_hotspots = hotspots(numpy_agg, kernel)
 
-    # cupy case
-    import cupy
-    cupy_agg = xr.DataArray(cupy.asarray(data))
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+def test_hotspot_gpu(data_hotspots):
+    data, kernel, expected_result = data_hotspots
+    cupy_agg = create_test_raster(data, backend='cupy')
     cupy_hotspots = hotspots(cupy_agg, kernel)
-
-    np.testing.assert_allclose(
-        numpy_hotspots.data, cupy_hotspots.data.get(), equal_nan=True)
+    general_output_checks(cupy_agg, cupy_hotspots, expected_result)
 
     # dask + cupy case not implemented
-    dask_cupy_agg = xr.DataArray(
-        da.from_array(cupy.asarray(data), chunks=(3, 3))
-    )
+    dask_cupy_agg = create_test_raster(data, backend='dask+cupy')
     with pytest.raises(NotImplementedError) as e_info:
         hotspots(dask_cupy_agg, kernel)
         assert e_info
diff --git a/xrspatial/tests/test_hillshade.py b/xrspatial/tests/test_hillshade.py
index bad0d66..aed0e56 100755
--- a/xrspatial/tests/test_hillshade.py
+++ b/xrspatial/tests/test_hillshade.py
@@ -3,39 +3,18 @@ import xarray as xr
 import numpy as np
 from numpy.testing import assert_allclose, assert_array_less
 
-import dask.array as da
-
 from xrspatial import hillshade
 from xrspatial.utils import doesnt_have_cuda
 from ..gpu_rtx import has_rtx
 
 from xrspatial.tests.general_checks import general_output_checks
+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
 
 
-elevation = np.asarray(
-    [[1432.6542, 1432.4764, 1432.4764, 1432.1207, 1431.9429, np.nan],
-     [1432.6542, 1432.6542, 1432.4764, 1432.2986, 1432.1207, np.nan],
-     [1432.832, 1432.6542, 1432.4764, 1432.2986, 1432.1207, np.nan],
-     [1432.832, 1432.6542, 1432.4764, 1432.4764, 1432.1207, np.nan],
-     [1432.832, 1432.6542, 1432.6542, 1432.4764, 1432.2986, np.nan],
-     [1432.832, 1432.6542, 1432.6542, 1432.4764, 1432.2986, np.nan],
-     [1432.832, 1432.832, 1432.6542, 1432.4764, 1432.4764, np.nan]],
-    dtype=np.float32)
-
-
-def _do_sparse_array(data_array):
-    import random
-    indx = list(zip(*np.where(data_array)))
-    pos = random.sample(range(data_array.size), data_array.size//2)
-    indx = np.asarray(indx)[pos]
-    r = indx[:, 0]
-    c = indx[:, 1]
-    data_half = data_array.copy()
-    data_half[r, c] = 0
-    return data_half
-
-
-def _do_gaussian_array():
+@pytest.fixture
+def data_gaussian():
     _x = np.linspace(0, 50, 101)
     _y = _x.copy()
     _mean = 25
@@ -47,12 +26,7 @@ def _do_gaussian_array():
     return gaussian
 
 
-data_random = np.random.random_sample((100, 100))
-data_random_sparse = _do_sparse_array(data_random)
-data_gaussian = _do_gaussian_array()
-
-
-def test_hillshade():
+def test_hillshade(data_gaussian):
     """
     Assert Simple Hillshade transfer function
     """
@@ -64,41 +38,27 @@ def test_hillshade():
     assert da_gaussian_shade[60, 60] > 0
 
 
-def test_hillshade_cpu():
-    attrs = {'res': (10.0, 10.0)}
-
-    small_numpy_based_data_array = xr.DataArray(elevation, attrs=attrs)
-    numpy_result = hillshade(small_numpy_based_data_array, name='numpy')
-
-    dask_data = da.from_array(elevation, chunks=(3, 3))
-    small_das_based_data_array = xr.DataArray(dask_data, attrs=attrs)
-    dask_result = hillshade(small_das_based_data_array, name='dask')
-
-    general_output_checks(small_das_based_data_array, dask_result)
-    np.testing.assert_allclose(
-        numpy_result.data, dask_result.data.compute(), equal_nan=True)
+@pytest.mark.parametrize("size", [(2, 4), (10, 15)])
+@pytest.mark.parametrize(
+    "dtype", [np.int32, np.int64, np.float32, np.float64])
+def test_hillshade_numpy_equals_dask_numpy(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, hillshade)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_hillshade_gpu_equals_cpu():
-
-    import cupy
-
-    small_da = xr.DataArray(elevation, attrs={'res': (10.0, 10.0)})
-    cpu = hillshade(small_da, name='numpy_result')
-
-    small_da_cupy = xr.DataArray(cupy.asarray(elevation),
-                                 attrs={'res': (10.0, 10.0)})
-    gpu = hillshade(small_da_cupy, name='cupy_result')
-
-    general_output_checks(small_da_cupy, gpu)
-    assert isinstance(gpu.data, cupy.ndarray)
-
-    assert np.isclose(cpu.data, gpu.data.get(), equal_nan=True).all()
+@pytest.mark.parametrize("size", [(2, 4), (10, 15)])
+@pytest.mark.parametrize(
+    "dtype", [np.int32, np.int64, np.float32, np.float64])
+def test_hillshade_gpu_equals_cpu(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, hillshade, rtol=1e-6)
 
 
 @pytest.mark.skipif(not has_rtx(), reason="RTX not available")
-def test_hillshade_rtx_with_shadows():
+def test_hillshade_rtx_with_shadows(data_gaussian):
     import cupy
 
     tall_gaussian = 400*data_gaussian
diff --git a/xrspatial/tests/test_multispectral.py b/xrspatial/tests/test_multispectral.py
index 52e4c6e..c5ac205 100755
--- a/xrspatial/tests/test_multispectral.py
+++ b/xrspatial/tests/test_multispectral.py
@@ -2,9 +2,6 @@ import pytest
 import xarray as xr
 import numpy as np
 
-import dask.array as da
-
-from xrspatial.utils import has_cuda
 from xrspatial.utils import doesnt_have_cuda
 
 from xrspatial.multispectral import arvi
@@ -20,27 +17,40 @@ from xrspatial.multispectral import sipi
 from xrspatial.multispectral import true_color
 
 from xrspatial.tests.general_checks import general_output_checks
-
-
-blue_data = np.array([[0, 9167, 9198, 9589.],
-                      [9154, 9122, 9202, 9475.],
-                      [9029, 9052, 9145, 9385.],
-                      [9153, 9134, 9157, 9389.],
-                      [9254, 9178, 9174, 9453.],
-                      [9348, 9222, 9201, 9413.],
-                      [9393, 9278, 9251, 9347.],
-                      [9486, 9293, np.nan, 9317.]])
-
-green_data = np.array([[0, 9929, 10056, 10620.],
-                       [9913, 9904, 10061, 10621.],
-                       [9853, 9874, 10116, 10429.],
-                       [9901, 9970, 10121, 10395.],
-                       [9954, 9945, 10068, 10512.],
-                       [9985, 9942, 10027, 10541.],
-                       [np.nan, 9980, 10062, 10371.],
-                       [10101, 9971, 10044, 10275.]])
-
-red_data = np.array([[0, 10301, 10454, 11325.],
+from xrspatial.tests.general_checks import create_test_raster
+
+
+@pytest.fixture
+def blue_data(backend):
+    data = np.array([[0, 9167, 9198, 9589.],
+                     [9154, 9122, 9202, 9475.],
+                     [9029, 9052, 9145, 9385.],
+                     [9153, 9134, 9157, 9389.],
+                     [9254, 9178, 9174, 9453.],
+                     [9348, 9222, 9201, 9413.],
+                     [9393, 9278, 9251, 9347.],
+                     [9486, 9293, np.nan, 9317.]])
+    agg = create_test_raster(data, backend=backend)
+    return agg
+
+
+@pytest.fixture
+def green_data(backend):
+    data = np.array([[0, 9929, 10056, 10620.],
+                     [9913, 9904, 10061, 10621.],
+                     [9853, 9874, 10116, 10429.],
+                     [9901, 9970, 10121, 10395.],
+                     [9954, 9945, 10068, 10512.],
+                     [9985, 9942, 10027, 10541.],
+                     [np.nan, 9980, 10062, 10371.],
+                     [10101, 9971, 10044, 10275.]])
+    agg = create_test_raster(data, backend=backend)
+    return agg
+
+
+@pytest.fixture
+def red_data(backend):
+    data = np.array([[0, 10301, 10454, 11325.],
                      [10353, 10269, 10501, np.nan],
                      [10184, 10106, 10445, 10961.],
                      [10349, 10230, 10299, 10844.],
@@ -48,8 +58,13 @@ red_data = np.array([[0, 10301, 10454, 11325.],
                      [10479, 10340, 10381, 11076.],
                      [10498, 10402, 10495, 10877.],
                      [10603, 10383, 10433, 10751.]])
+    agg = create_test_raster(data, backend=backend)
+    return agg
 
-nir_data = np.array([[0, 15928, 16135, 16411.],
+
+@pytest.fixture
+def nir_data(backend):
+    data = np.array([[0, 15928, 16135, 16411.],
                      [15588, 15881, 16253, 16651.],
                      [16175, 16486, 17038, 17084.],
                      [15671, 16596, 17511, 17525.],
@@ -57,8 +72,13 @@ nir_data = np.array([[0, 15928, 16135, 16411.],
                      [15317, 15782, 16322, 17133.],
                      [15168, 15529, 16011, 16600.],
                      [15072, 15496, 15983, 16477.]])
+    agg = create_test_raster(data, backend=backend)
+    return agg
+
 
-tir_data = np.array([[0, 10512, 10517, 10527.],
+@pytest.fixture
+def tir_data(backend):
+    data = np.array([[0, 10512, 10517, 10527.],
                      [10511, 10504, 10502, 10504.],
                      [10522, 10507, 10497, 10491.],
                      [10543, 10514, 10498, 10486.],
@@ -66,170 +86,187 @@ tir_data = np.array([[0, 10512, 10517, 10527.],
                      [10592, 10558, 10527, 10504.],
                      [10629, 10598, 10567, 10536.],
                      [10664, 10639, 10612, 10587.]])
-
-swir1_data = np.array([[0, np.nan, 17194, 18163.],
-                       [16974, 16871, 17123, 18304.],
-                       [16680, 16437, 16474, 17519.],
-                       [17004, 16453, 16001, 16800.],
-                       [17230, 16906, 16442, 16840.],
-                       [17237, 16969, 16784, 17461.],
-                       [17417, 17079, 17173, 17679.],
-                       [17621, 17205, 17163, 17362.]])
-
-swir2_data = np.array([[0, 13175, 13558, 14952.],
-                       [13291, 13159, 13516, 15029.],
-                       [12924, 12676, np.nan, 14009.],
-                       [13294, 12728, 12370, 13289.],
-                       [13507, 13163, 12763, 13499.],
-                       [13570, 13219, 13048, 14145.],
-                       [13770, 13393, 13472, 14249.],
-                       [14148, 13489, 13483, 13893.]])
-
-arvi_expected_results = np.array([
-    [np.nan, 0.09832155, 0.0956943, 0.0688592],
-    [0.08880479, 0.09804352, 0.09585208, np.nan],
-    [0.10611779, 0.1164153, 0.11244237, 0.09396376],
-    [0.0906375, 0.11409396, 0.12842213, 0.10752644],
-    [0.08580945, 0.09740005, 0.1179347, 0.10302287],
-    [0.08125288, 0.09465021, 0.1028627, 0.09022958],
-    [0.07825362, 0.08776391, 0.09236357, 0.08790172],
-    [0.07324535, 0.08831083, np.nan, 0.09074763]], dtype=np.float32)
-
-evi_expected_results = np.array([
-    [0., 1.5661007, 1.4382279, 1.0217365],
-    [1.4458131, 1.544984, 1.4036115, np.nan],
-    [1.5662745, 1.7274992, 1.4820393, 1.2281862],
-    [1.4591216, 1.6802154, 1.6963824, 1.3721503],
-    [1.4635549, 1.5457553, 1.6425549, 1.3112202],
-    [1.4965355, 1.5713791, 1.5468937, 1.1654801],
-    [1.5143654, 1.5337442, 1.4365331, 1.2165724],
-    [1.4805857, 1.5785736, np.nan, 1.2888849]], dtype=np.float32)
-
-nbr_expected_results = np.array([
-    [np.nan, 0.09459506, 0.08678813, 0.04651979],
-    [0.07953876, 0.09373278, 0.09194128, 0.0511995],
-    [0.11172205, 0.13064948, np.nan, 0.09889686],
-    [0.08206456, 0.1319056, 0.17204913, 0.13746998],
-    [0.06941334, 0.09529537, 0.1424444, 0.13044319],
-    [0.06047703, 0.08837626, 0.11147429, 0.09553041],
-    [0.04831018, 0.07385381, 0.08611742, 0.07620993],
-    [0.03162218, 0.06924271, 0.08484355, 0.08508396]], dtype=np.float32)
-
-nbr2_expected_results = np.array([
-    [np.nan, np.nan, 0.11823621, 0.09696512],
-    [0.12169173, 0.12360972, 0.11772577, 0.09825099],
-    [0.12687474, 0.12918627, np.nan, 0.11132962],
-    [0.12245033, 0.12765156, 0.1279828, 0.11668716],
-    [0.12112438, 0.12448036, 0.12597159, 0.11012229],
-    [0.11903139, 0.12422155, 0.12523465, 0.10491679],
-    [0.11693975, 0.12096351, 0.12077011, 0.10742921],
-    [0.10932041, 0.121066, 0.12008093, 0.11099024]], dtype=np.float32)
-
-ndvi_expected_results = np.array([
-    [np.nan, 0.21453354, 0.21365978, 0.1833718],
-    [0.20180409, 0.21460803, 0.21499589, np.nan],
-    [0.2272848, 0.23992178, 0.23989375, 0.21832769],
-    [0.20453498, 0.23730709, 0.25933117, 0.23550354],
-    [0.19620839, 0.21352422, 0.24427369, 0.23030005],
-    [0.18754846, 0.20833014, 0.22248437, 0.2147187],
-    [0.18195277, 0.19771701, 0.20810382, 0.20828329],
-    [0.17406037, 0.19757332, 0.21009994, 0.21029823]], dtype=np.float32)
-
-ndmi_expected_results = np.array([
-    [np.nan, np.nan, -0.03177413, -0.05067392],
-    [-0.04256495, -0.03022716, -0.02606663, -0.04728937],
-    [-0.01537057,  0.00148832,  0.01682979, -0.01257116],
-    [-0.04079571,  0.00432691,  0.04505849,  0.02112163],
-    [-0.05214949, -0.02953535,  0.01677381,  0.02061706],
-    [-0.05897893, -0.03624317, -0.01395517, -0.00948141],
-    [-0.06901949, -0.04753435, -0.03501688, -0.031477],
-    [-0.07796776, -0.0522614, -0.03560007, -0.02615326]], dtype=np.float32)
-
-savi_expected_results = np.array([
-    [0., 0.10726268, 0.10682587, 0.09168259],
-    [0.10089815, 0.10729991, 0.10749393, np.nan],
-    [0.11363809, 0.11995638, 0.11994251, 0.10915995],
-    [0.10226355, 0.11864913, 0.12966092, 0.11774762],
-    [0.09810041, 0.10675804, 0.12213238, 0.11514599],
-    [0.09377059, 0.10416108, 0.11123802, 0.10735555],
-    [0.09097284, 0.0988547, 0.10404798, 0.10413785],
-    [0.0870268, 0.09878284, 0.105046, 0.10514525]], dtype=np.float32)
-
-gci_expected_results = np.array([
-    [np.nan, 0.60418975, 0.6045147, 0.5452919],
-    [0.57248056, 0.6034935, 0.6154458, 0.5677431],
-    [0.64163196, 0.66963744, 0.6842626, 0.63812447],
-    [0.5827694, 0.66459376, 0.730165, 0.6859067],
-    [0.55937314, 0.6024133, 0.6888161, 0.6694254],
-    [0.534001, 0.58740693, 0.62780493, 0.62536764],
-    [np.nan, 0.55601203, 0.5912343, 0.6006171],
-    [0.4921295, 0.5541069, 0.5912983, 0.603601]], dtype=np.float32)
-
-sipi_expected_results = np.array([
-    [np.nan, 1.2015283, 1.2210878, 1.3413291],
-    [1.2290354, 1.2043835, 1.2258345, np.nan],
-    [1.1927892, 1.1652038, 1.1971788, 1.2573901],
-    [1.2247275, 1.1721647, 1.1583472, 1.2177818],
-    [1.2309505, 1.2050642, 1.1727082, 1.2322679],
-    [1.2337743, 1.2054392, 1.1986197, 1.2745583],
-    [1.2366167, 1.2192315, 1.2255257, 1.2673423],
-    [1.2499441, 1.2131821, np.nan, 1.2504367]], dtype=np.float32)
-
-ebbi_expected_results = np.array([
-    [np.nan, np.nan, 4.0488696, 4.0370474],
-    [3.9937027, 3.9902349, 3.9841716, np.nan],
-    [3.9386337, 3.8569257, 3.6711047, 3.918455],
-    [4.0096908, 3.7895138, 3.5027769, 3.6056597],
-    [4.0786624, np.nan, 3.724852, 3.5452912],
-    [4.0510664, 3.9954765, 3.8744915, 3.8181543],
-    [4.131501, 4.013487, 4.009527, 4.049455],
-    [4.172874, 4.08833, 4.038202, 3.954431]], dtype=np.float32)
-
-
-def _do_gaussian_array():
-    _x = np.linspace(0, 50, 101)
-    _y = _x.copy()
-    _mean = 25
-    _sdev = 5
-    X, Y = np.meshgrid(_x, _y, sparse=True)
-    x_fac = -np.power(X-_mean, 2)
-    y_fac = -np.power(Y-_mean, 2)
-    gaussian = np.exp((x_fac+y_fac)/(2*_sdev**2)) / (2.5*_sdev)
-    return gaussian
-
-
-data_gaussian = _do_gaussian_array()
-
-
-def create_test_arr(arr, backend='numpy'):
-
-    y, x = arr.shape
-    raster = xr.DataArray(arr, dims=['y', 'x'])
-
-    if backend == 'numpy':
-        raster['y'] = np.linspace(0, y, y)
-        raster['x'] = np.linspace(0, x, x)
-        return raster
-
-    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=(3, 3))
-
-    return raster
+    agg = create_test_raster(data, backend=backend)
+    return agg
+
+
+@pytest.fixture
+def swir1_data(backend):
+    data = np.array([[0, np.nan, 17194, 18163.],
+                     [16974, 16871, 17123, 18304.],
+                     [16680, 16437, 16474, 17519.],
+                     [17004, 16453, 16001, 16800.],
+                     [17230, 16906, 16442, 16840.],
+                     [17237, 16969, 16784, 17461.],
+                     [17417, 17079, 17173, 17679.],
+                     [17621, 17205, 17163, 17362.]])
+    agg = create_test_raster(data, backend=backend)
+    return agg
+
+
+@pytest.fixture
+def swir2_data(backend):
+    data = np.array([[0, 13175, 13558, 14952.],
+                     [13291, 13159, 13516, 15029.],
+                     [12924, 12676, np.nan, 14009.],
+                     [13294, 12728, 12370, 13289.],
+                     [13507, 13163, 12763, 13499.],
+                     [13570, 13219, 13048, 14145.],
+                     [13770, 13393, 13472, 14249.],
+                     [14148, 13489, 13483, 13893.]])
+    agg = create_test_raster(data, backend=backend)
+    return agg
+
+
+@pytest.fixture
+def result_arvi():
+    result = np.array([
+        [np.nan, 0.09832155, 0.0956943, 0.0688592],
+        [0.08880479, 0.09804352, 0.09585208, np.nan],
+        [0.10611779, 0.1164153, 0.11244237, 0.09396376],
+        [0.0906375, 0.11409396, 0.12842213, 0.10752644],
+        [0.08580945, 0.09740005, 0.1179347, 0.10302287],
+        [0.08125288, 0.09465021, 0.1028627, 0.09022958],
+        [0.07825362, 0.08776391, 0.09236357, 0.08790172],
+        [0.07324535, 0.08831083, np.nan, 0.09074763]], dtype=np.float32)
+    return result
+
+
+@pytest.fixture
+def result_evi():
+    result = np.array([
+        [0., 1.5661007, 1.4382279, 1.0217365],
+        [1.4458131, 1.544984, 1.4036115, np.nan],
+        [1.5662745, 1.7274992, 1.4820393, 1.2281862],
+        [1.4591216, 1.6802154, 1.6963824, 1.3721503],
+        [1.4635549, 1.5457553, 1.6425549, 1.3112202],
+        [1.4965355, 1.5713791, 1.5468937, 1.1654801],
+        [1.5143654, 1.5337442, 1.4365331, 1.2165724],
+        [1.4805857, 1.5785736, np.nan, 1.2888849]], dtype=np.float32)
+    return result
+
+
+@pytest.fixture
+def result_nbr():
+    result = np.array([
+        [np.nan, 0.09459506, 0.08678813, 0.04651979],
+        [0.07953876, 0.09373278, 0.09194128, 0.0511995],
+        [0.11172205, 0.13064948, np.nan, 0.09889686],
+        [0.08206456, 0.1319056, 0.17204913, 0.13746998],
+        [0.06941334, 0.09529537, 0.1424444, 0.13044319],
+        [0.06047703, 0.08837626, 0.11147429, 0.09553041],
+        [0.04831018, 0.07385381, 0.08611742, 0.07620993],
+        [0.03162218, 0.06924271, 0.08484355, 0.08508396]], dtype=np.float32)
+    return result
+
+
+@pytest.fixture
+def result_nbr2():
+    result = np.array([
+        [np.nan, np.nan, 0.11823621, 0.09696512],
+        [0.12169173, 0.12360972, 0.11772577, 0.09825099],
+        [0.12687474, 0.12918627, np.nan, 0.11132962],
+        [0.12245033, 0.12765156, 0.1279828, 0.11668716],
+        [0.12112438, 0.12448036, 0.12597159, 0.11012229],
+        [0.11903139, 0.12422155, 0.12523465, 0.10491679],
+        [0.11693975, 0.12096351, 0.12077011, 0.10742921],
+        [0.10932041, 0.121066, 0.12008093, 0.11099024]], dtype=np.float32)
+    return result
+
+
+@pytest.fixture
+def result_ndvi():
+    result = np.array([
+        [np.nan, 0.21453354, 0.21365978, 0.1833718],
+        [0.20180409, 0.21460803, 0.21499589, np.nan],
+        [0.2272848, 0.23992178, 0.23989375, 0.21832769],
+        [0.20453498, 0.23730709, 0.25933117, 0.23550354],
+        [0.19620839, 0.21352422, 0.24427369, 0.23030005],
+        [0.18754846, 0.20833014, 0.22248437, 0.2147187],
+        [0.18195277, 0.19771701, 0.20810382, 0.20828329],
+        [0.17406037, 0.19757332, 0.21009994, 0.21029823]], dtype=np.float32)
+    return result
+
+
+@pytest.fixture
+def result_ndmi():
+    result = np.array([
+        [np.nan, np.nan, -0.03177413, -0.05067392],
+        [-0.04256495, -0.03022716, -0.02606663, -0.04728937],
+        [-0.01537057,  0.00148832,  0.01682979, -0.01257116],
+        [-0.04079571,  0.00432691,  0.04505849,  0.02112163],
+        [-0.05214949, -0.02953535,  0.01677381,  0.02061706],
+        [-0.05897893, -0.03624317, -0.01395517, -0.00948141],
+        [-0.06901949, -0.04753435, -0.03501688, -0.031477],
+        [-0.07796776, -0.0522614, -0.03560007, -0.02615326]], dtype=np.float32)
+    return result
+
+
+@pytest.fixture
+def result_savi():
+    result = np.array([
+        [0., 0.10726268, 0.10682587, 0.09168259],
+        [0.10089815, 0.10729991, 0.10749393, np.nan],
+        [0.11363809, 0.11995638, 0.11994251, 0.10915995],
+        [0.10226355, 0.11864913, 0.12966092, 0.11774762],
+        [0.09810041, 0.10675804, 0.12213238, 0.11514599],
+        [0.09377059, 0.10416108, 0.11123802, 0.10735555],
+        [0.09097284, 0.0988547, 0.10404798, 0.10413785],
+        [0.0870268, 0.09878284, 0.105046, 0.10514525]], dtype=np.float32)
+    return result
+
+
+@pytest.fixture
+def result_gci():
+    result = np.array([
+        [np.nan, 0.60418975, 0.6045147, 0.5452919],
+        [0.57248056, 0.6034935, 0.6154458, 0.5677431],
+        [0.64163196, 0.66963744, 0.6842626, 0.63812447],
+        [0.5827694, 0.66459376, 0.730165, 0.6859067],
+        [0.55937314, 0.6024133, 0.6888161, 0.6694254],
+        [0.534001, 0.58740693, 0.62780493, 0.62536764],
+        [np.nan, 0.55601203, 0.5912343, 0.6006171],
+        [0.4921295, 0.5541069, 0.5912983, 0.603601]], dtype=np.float32)
+    return result
+
+
+@pytest.fixture
+def result_sipi():
+    result = np.array([
+        [np.nan, 1.2015283, 1.2210878, 1.3413291],
+        [1.2290354, 1.2043835, 1.2258345, np.nan],
+        [1.1927892, 1.1652038, 1.1971788, 1.2573901],
+        [1.2247275, 1.1721647, 1.1583472, 1.2177818],
+        [1.2309505, 1.2050642, 1.1727082, 1.2322679],
+        [1.2337743, 1.2054392, 1.1986197, 1.2745583],
+        [1.2366167, 1.2192315, 1.2255257, 1.2673423],
+        [1.2499441, 1.2131821, np.nan, 1.2504367]], dtype=np.float32)
+    return result
+
+
+@pytest.fixture
+def result_ebbi():
+    result = np.array([
+        [np.nan, np.nan, 4.0488696, 4.0370474],
+        [3.9937027, 3.9902349, 3.9841716, np.nan],
+        [3.9386337, 3.8569257, 3.6711047, 3.918455],
+        [4.0096908, 3.7895138, 3.5027769, 3.6056597],
+        [4.0786624, np.nan, 3.724852, 3.5452912],
+        [4.0510664, 3.9954765, 3.8744915, 3.8181543],
+        [4.131501, 4.013487, 4.009527, 4.049455],
+        [4.172874, 4.08833, 4.038202, 3.954431]], dtype=np.float32)
+    return result
 
 
 # NDVI -------------
-def test_ndvi_numpy_contains_valid_values():
+def test_ndvi_data_contains_valid_values():
     _x = np.mgrid[1:0:21j]
     a, b = np.meshgrid(_x, _x)
-    red_numpy = a*b
-    nir_numpy = (a*b)[::-1, ::-1]
+    red_data = a*b
+    nir_data = (a*b)[::-1, ::-1]
 
-    da_nir = xr.DataArray(nir_numpy, dims=['y', 'x'])
-    da_red = xr.DataArray(red_numpy, dims=['y', 'x'])
+    da_nir = xr.DataArray(nir_data, dims=['y', 'x'])
+    da_red = xr.DataArray(red_data, dims=['y', 'x'])
 
     da_ndvi = ndvi(da_nir, da_red)
 
@@ -244,337 +281,170 @@ def test_ndvi_numpy_contains_valid_values():
     assert da_ndvi[15, 10] == da_ndvi[10, 15] == 0.5
 
 
-def test_ndvi_cpu():
-
-    # vanilla numpy version
-    nir_numpy = create_test_arr(nir_data)
-    red_numpy = create_test_arr(red_data)
-    numpy_result = ndvi(nir_numpy, red_numpy)
-    general_output_checks(nir_numpy, numpy_result, ndvi_expected_results)
-
-    # dask
-    nir_dask = create_test_arr(nir_data, backend='dask')
-    red_dask = create_test_arr(red_data, backend='dask')
-    dask_result = ndvi(nir_dask, red_dask)
-    general_output_checks(nir_dask, dask_result, ndvi_expected_results)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_ndvi_cpu(nir_data, red_data, result_ndvi):
+    result = ndvi(nir_data, red_data)
+    general_output_checks(nir_data, result, result_ndvi, verify_dtype=True)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_ndvi_gpu():
-    # cupy
-    nir_cupy = create_test_arr(nir_data, backend='cupy')
-    red_cupy = create_test_arr(red_data, backend='cupy')
-    cupy_result = ndvi(nir_cupy, red_cupy)
-    general_output_checks(nir_cupy, cupy_result, ndvi_expected_results)
-
-    # dask + cupy
-    nir_dask = create_test_arr(nir_data, backend='dask+cupy')
-    red_dask = create_test_arr(red_data, backend='dask+cupy')
-    dask_cupy_result = ndvi(nir_dask, red_dask)
-    general_output_checks(nir_dask, dask_cupy_result, ndvi_expected_results)
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_ndvi_gpu(nir_data, red_data, result_ndvi):
+    result = ndvi(nir_data, red_data)
+    general_output_checks(nir_data, result, result_ndvi, verify_dtype=True)
 
 
 # SAVI -------------
-def test_savi_cpu():
-    nir_numpy = create_test_arr(nir_data)
-    red_numpy = create_test_arr(red_data)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_savi_zero_soil_factor_cpu(nir_data, red_data, result_ndvi):
+    # savi should be same as ndvi at soil_factor=0
+    result_savi = savi(nir_data, red_data, soil_factor=0.0)
+    general_output_checks(nir_data, result_savi, result_ndvi, verify_dtype=True)
 
+
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_savi_zero_soil_factor_gpu(nir_data, red_data, result_ndvi):
     # savi should be same as ndvi at soil_factor=0
-    result_savi = savi(nir_numpy, red_numpy, soil_factor=0.0)
-    result_ndvi = ndvi(nir_numpy, red_numpy)
-    assert np.isclose(result_savi.data, result_ndvi.data, equal_nan=True).all()
+    result_savi = savi(nir_data, red_data, soil_factor=0.0)
+    general_output_checks(nir_data, result_savi, result_ndvi, verify_dtype=True)
 
-    # test default savi where soil_factor = 1.0
-    numpy_result = savi(nir_numpy, red_numpy, soil_factor=1.0)
-    general_output_checks(nir_numpy, numpy_result, savi_expected_results)
 
-    # dask
-    nir_dask = create_test_arr(nir_data, backend='dask')
-    red_dask = create_test_arr(red_data, backend='dask')
-    dask_result = savi(nir_dask, red_dask)
-    general_output_checks(nir_dask, dask_result, savi_expected_results)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_savi_cpu(nir_data, red_data, result_savi):
+    # test default savi where soil_factor = 1.0
+    result = savi(nir_data, red_data, soil_factor=1.0)
+    general_output_checks(nir_data, result, result_savi)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_savi_gpu():
-    # cupy
-    nir_cupy = create_test_arr(nir_data, backend='cupy')
-    red_cupy = create_test_arr(red_data, backend='cupy')
-    cupy_result = savi(nir_cupy, red_cupy)
-    general_output_checks(nir_cupy, cupy_result, savi_expected_results)
-
-    # dask + cupy
-    nir_dask_cupy = create_test_arr(nir_data, backend='dask+cupy')
-    red_dask_cupy = create_test_arr(red_data, backend='dask+cupy')
-    dask_cupy_result = savi(nir_dask_cupy, red_dask_cupy)
-    general_output_checks(
-        nir_dask_cupy, dask_cupy_result, savi_expected_results)
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_savi_gpu(nir_data, red_data, result_savi):
+    # test default savi where soil_factor = 1.0
+    result = savi(nir_data, red_data, soil_factor=1.0)
+    general_output_checks(nir_data, result, result_savi)
 
 
 # arvi -------------
-def test_arvi_cpu():
-    nir_numpy = create_test_arr(nir_data)
-    red_numpy = create_test_arr(red_data)
-    blue_numpy = create_test_arr(blue_data)
-    numpy_result = arvi(nir_numpy, red_numpy, blue_numpy)
-    general_output_checks(nir_numpy, numpy_result, arvi_expected_results)
-
-    # dask
-    nir_dask = create_test_arr(nir_data, backend='dask')
-    red_dask = create_test_arr(red_data, backend='dask')
-    blue_dask = create_test_arr(blue_data, backend='dask')
-    dask_result = arvi(nir_dask, red_dask, blue_dask)
-    general_output_checks(nir_dask, dask_result, arvi_expected_results)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_arvi_cpu(nir_data, red_data, blue_data, result_arvi):
+    result = arvi(nir_data, red_data, blue_data)
+    general_output_checks(nir_data, result, result_arvi)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_arvi_gpu():
-    # cupy
-    nir_cupy = create_test_arr(nir_data, backend='cupy')
-    red_cupy = create_test_arr(red_data, backend='cupy')
-    blue_cupy = create_test_arr(blue_data, backend='cupy')
-    cupy_result = arvi(nir_cupy, red_cupy, blue_cupy)
-    general_output_checks(nir_cupy, cupy_result, arvi_expected_results)
-
-    # dask + cupy
-    nir_dask_cupy = create_test_arr(nir_data, backend='dask+cupy')
-    red_dask_cupy = create_test_arr(red_data, backend='dask+cupy')
-    blue_dask_cupy = create_test_arr(blue_data, backend='dask+cupy')
-    dask_cupy_result = arvi(nir_dask_cupy, red_dask_cupy, blue_dask_cupy)
-    general_output_checks(
-        nir_dask_cupy, dask_cupy_result, arvi_expected_results
-    )
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_arvi_gpu(nir_data, red_data, blue_data, result_arvi):
+    result = arvi(nir_data, red_data, blue_data)
+    general_output_checks(nir_data, result, result_arvi)
 
 
 # EVI -------------
-def test_evi_cpu():
-    nir_numpy = create_test_arr(nir_data)
-    red_numpy = create_test_arr(red_data)
-    blue_numpy = create_test_arr(blue_data)
-    numpy_result = evi(nir_numpy, red_numpy, blue_numpy)
-    general_output_checks(nir_numpy, numpy_result, evi_expected_results)
-
-    # dask
-    nir_dask = create_test_arr(nir_data, backend='dask')
-    red_dask = create_test_arr(red_data, backend='dask')
-    blue_dask = create_test_arr(blue_data, backend='dask')
-    dask_result = evi(nir_dask, red_dask, blue_dask)
-    general_output_checks(nir_dask, dask_result, evi_expected_results)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_evi_cpu(nir_data, red_data, blue_data, result_evi):
+    result = evi(nir_data, red_data, blue_data)
+    general_output_checks(nir_data, result, result_evi)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_evi_gpu():
-    # cupy
-    nir_cupy = create_test_arr(nir_data, backend='cupy')
-    red_cupy = create_test_arr(red_data, backend='cupy')
-    blue_cupy = create_test_arr(blue_data, backend='cupy')
-    cupy_result = evi(nir_cupy, red_cupy, blue_cupy)
-    general_output_checks(nir_cupy, cupy_result, evi_expected_results)
-
-    # dask + cupy
-    nir_dask_cupy = create_test_arr(nir_data, backend='dask+cupy')
-    red_dask_cupy = create_test_arr(red_data, backend='dask+cupy')
-    blue_dask_cupy = create_test_arr(blue_data, backend='dask+cupy')
-    dask_cupy_result = evi(nir_dask_cupy, red_dask_cupy, blue_dask_cupy)
-    general_output_checks(
-        nir_dask_cupy, dask_cupy_result, evi_expected_results)
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_evi_gpu(nir_data, red_data, blue_data, result_evi):
+    result = evi(nir_data, red_data, blue_data)
+    general_output_checks(nir_data, result, result_evi)
 
 
 # GCI -------------
-def test_gci_cpu():
-    # vanilla numpy version
-    nir_numpy = create_test_arr(nir_data)
-    green_numpy = create_test_arr(green_data)
-    numpy_result = gci(nir_numpy, green_numpy)
-    general_output_checks(nir_numpy, numpy_result, gci_expected_results)
-
-    # dask
-    nir_dask = create_test_arr(nir_data, backend='dask')
-    green_dask = create_test_arr(green_data, backend='dask')
-    dask_result = gci(nir_dask, green_dask)
-    general_output_checks(nir_dask, dask_result, gci_expected_results)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_gci_cpu(nir_data, green_data, result_gci):
+    result = gci(nir_data, green_data)
+    general_output_checks(nir_data, result, result_gci)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_gci_gpu():
-    # cupy
-    nir_cupy = create_test_arr(nir_data, backend='cupy')
-    green_cupy = create_test_arr(green_data, backend='cupy')
-    cupy_result = gci(nir_cupy, green_cupy)
-    general_output_checks(nir_cupy, cupy_result, gci_expected_results)
-
-    # dask + cupy
-    nir_dask_cupy = create_test_arr(nir_data, backend='dask+cupy')
-    green_dask_cupy = create_test_arr(green_data, backend='dask+cupy')
-    dask_cupy_result = gci(nir_dask_cupy, green_dask_cupy)
-    general_output_checks(
-        nir_dask_cupy, dask_cupy_result, gci_expected_results)
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_gci_gpu(nir_data, green_data, result_gci):
+    result = gci(nir_data, green_data)
+    general_output_checks(nir_data, result, result_gci)
 
 
 # SIPI -------------
-def test_sipi_cpu():
-    nir_numpy = create_test_arr(nir_data)
-    red_numpy = create_test_arr(red_data)
-    blue_numpy = create_test_arr(blue_data)
-    numpy_result = sipi(nir_numpy, red_numpy, blue_numpy)
-    general_output_checks(nir_numpy, numpy_result, sipi_expected_results)
-
-    # dask
-    nir_dask = create_test_arr(nir_data, backend='dask')
-    red_dask = create_test_arr(red_data, backend='dask')
-    blue_dask = create_test_arr(blue_data, backend='dask')
-    dask_result = sipi(nir_dask, red_dask, blue_dask)
-    general_output_checks(nir_dask, dask_result, sipi_expected_results)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_sipi_cpu(nir_data, red_data, blue_data, result_sipi):
+    result = sipi(nir_data, red_data, blue_data)
+    general_output_checks(nir_data, result, result_sipi)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_sipi_gpu():
-    # cupy
-    nir_cupy = create_test_arr(nir_data, backend='cupy')
-    red_cupy = create_test_arr(red_data, backend='cupy')
-    blue_cupy = create_test_arr(blue_data, backend='cupy')
-    cupy_result = sipi(nir_cupy, red_cupy, blue_cupy)
-    general_output_checks(nir_cupy, cupy_result, sipi_expected_results)
-
-    # dask + cupy
-    nir_dask_cupy = create_test_arr(nir_data, backend='dask+cupy')
-    red_dask_cupy = create_test_arr(red_data, backend='dask+cupy')
-    blue_dask_cupy = create_test_arr(blue_data, backend='dask+cupy')
-    dask_cupy_result = sipi(nir_dask_cupy, red_dask_cupy, blue_dask_cupy)
-    general_output_checks(
-        nir_dask_cupy, dask_cupy_result, sipi_expected_results)
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_sipi_gpu(nir_data, red_data, blue_data, result_sipi):
+    result = sipi(nir_data, red_data, blue_data)
+    general_output_checks(nir_data, result, result_sipi)
 
 
 # NBR -------------
-def test_nbr_cpu():
-    nir_numpy = create_test_arr(nir_data)
-    swir_numpy = create_test_arr(swir2_data)
-    numpy_result = nbr(nir_numpy, swir_numpy)
-    general_output_checks(nir_numpy, numpy_result, nbr_expected_results)
-
-    # dask
-    nir_dask = create_test_arr(nir_data, backend='dask')
-    swir_dask = create_test_arr(swir2_data, backend='dask')
-    dask_result = nbr(nir_dask, swir_dask)
-    general_output_checks(nir_dask, dask_result, nbr_expected_results)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_nbr_cpu(nir_data, swir2_data, result_nbr):
+    result = nbr(nir_data, swir2_data)
+    general_output_checks(nir_data, result, result_nbr)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_nbr_gpu():
-    # cupy
-    nir_cupy = create_test_arr(nir_data, backend='cupy')
-    swir_cupy = create_test_arr(swir2_data, backend='cupy')
-    cupy_result = nbr(nir_cupy, swir_cupy)
-    general_output_checks(nir_cupy, cupy_result, nbr_expected_results)
-
-    # dask + cupy
-    nir_dask_cupy = create_test_arr(nir_data, backend='dask+cupy')
-    swir_dask_cupy = create_test_arr(swir2_data, backend='dask+cupy')
-    dask_cupy_result = nbr(nir_dask_cupy, swir_dask_cupy)
-    general_output_checks(
-        nir_dask_cupy, dask_cupy_result, nbr_expected_results)
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_nbr_gpu(nir_data, swir2_data, result_nbr):
+    result = nbr(nir_data, swir2_data)
+    general_output_checks(nir_data, result, result_nbr)
 
 
 # NBR2 -------------
-def test_nbr2_cpu():
-    swir1_numpy = create_test_arr(swir1_data)
-    swir2_numpy = create_test_arr(swir2_data)
-    numpy_result = nbr2(swir1_numpy, swir2_numpy)
-    general_output_checks(swir1_numpy, numpy_result, nbr2_expected_results)
-
-    # dask
-    swir1_dask = create_test_arr(swir1_data, backend='dask')
-    swir2_dask = create_test_arr(swir2_data, backend='dask')
-    dask_result = nbr2(swir1_dask, swir2_dask)
-    general_output_checks(swir1_dask, dask_result, nbr2_expected_results)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_nbr2_cpu(swir1_data, swir2_data, result_nbr2):
+    result = nbr2(swir1_data, swir2_data)
+    general_output_checks(swir1_data, result, result_nbr2)
 
 
-@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Dnbr2ce not Available")
-def test_nbr2_gpu():
-    # cupy
-    swir1_cupy = create_test_arr(swir1_data, backend='cupy')
-    swir2_cupy = create_test_arr(swir2_data, backend='cupy')
-    cupy_result = nbr2(swir1_cupy, swir2_cupy)
-    general_output_checks(swir2_cupy, cupy_result, nbr2_expected_results)
-
-    # dask + cupy
-    swir1_dask_cupy = create_test_arr(swir1_data, backend='dask+cupy')
-    swir2_dask_cupy = create_test_arr(swir2_data, backend='dask+cupy')
-    dask_cupy_result = nbr2(swir1_dask_cupy, swir2_dask_cupy)
-    general_output_checks(
-        swir1_dask_cupy, dask_cupy_result, nbr2_expected_results)
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_nbr2_gpu(swir1_data, swir2_data, result_nbr2):
+    result = nbr2(swir1_data, swir2_data)
+    general_output_checks(swir1_data, result, result_nbr2)
 
 
 # NDMI -------------
-def test_ndmi_cpu():
-    nir_numpy = create_test_arr(nir_data)
-    swir1_numpy = create_test_arr(swir1_data)
-    numpy_result = ndmi(nir_numpy, swir1_numpy)
-    general_output_checks(nir_numpy, numpy_result, ndmi_expected_results)
-
-    # dask
-    nir_dask = create_test_arr(nir_data, backend='dask')
-    swir1_dask = create_test_arr(swir1_data, backend='dask')
-    dask_result = ndmi(nir_dask, swir1_dask)
-    general_output_checks(nir_dask, dask_result, ndmi_expected_results)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_ndmi_cpu(nir_data, swir1_data, result_ndmi):
+    result = ndmi(nir_data, swir1_data)
+    general_output_checks(nir_data, result, result_ndmi)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_ndmi_gpu():
-    # cupy
-    nir_cupy = create_test_arr(nir_data, backend='cupy')
-    swir1_cupy = create_test_arr(swir1_data, backend='cupy')
-    cupy_result = ndmi(nir_cupy, swir1_cupy)
-    general_output_checks(nir_cupy, cupy_result, ndmi_expected_results)
-
-    # dask + cupy
-    nir_dask_cupy = create_test_arr(nir_data, backend='dask+cupy')
-    swir1_dask_cupy = create_test_arr(swir1_data, backend='dask+cupy')
-    dask_cupy_result = ndmi(nir_dask_cupy, swir1_dask_cupy)
-    general_output_checks(
-        nir_dask_cupy, dask_cupy_result, ndmi_expected_results)
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_ndmi_gpu(nir_data, swir1_data, result_ndmi):
+    result = ndmi(nir_data, swir1_data)
+    general_output_checks(nir_data, result, result_ndmi)
 
 
 # EBBI -------------
-def test_ebbi_cpu():
-    # vanilla numpy version
-    red_numpy = create_test_arr(red_data)
-    swir_numpy = create_test_arr(swir1_data)
-    tir_numpy = create_test_arr(tir_data)
-    numpy_result = ebbi(red_numpy, swir_numpy, tir_numpy)
-    general_output_checks(red_numpy, numpy_result, ebbi_expected_results)
-
-    # dask
-    red_dask = create_test_arr(red_data, backend='dask')
-    swir_dask = create_test_arr(swir1_data, backend='dask')
-    tir_dask = create_test_arr(tir_data, backend='dask')
-    dask_result = ebbi(red_dask, swir_dask, tir_dask)
-    general_output_checks(red_dask, dask_result, ebbi_expected_results)
+@pytest.mark.parametrize("backend", ["numpy", "dask+numpy"])
+def test_ebbi_cpu(red_data, swir1_data, tir_data, result_ebbi):
+    result = ebbi(red_data, swir1_data, tir_data)
+    general_output_checks(red_data, result, result_ebbi)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_ebbi_gpu():
-    # cupy
-    red_cupy = create_test_arr(red_data, backend='cupy')
-    swir_cupy = create_test_arr(swir1_data, backend='cupy')
-    tir_cupy = create_test_arr(tir_data, backend='cupy')
-    cupy_result = ebbi(red_cupy, swir_cupy, tir_cupy)
-    general_output_checks(red_cupy, cupy_result, ebbi_expected_results)
-
-    # dask + cupy
-    red_dask_cupy = create_test_arr(red_data, backend='dask+cupy')
-    swir_dask_cupy = create_test_arr(swir1_data, backend='dask+cupy')
-    tir_dask_cupy = create_test_arr(tir_data, backend='dask+cupy')
-    dask_cupy_result = ebbi(red_dask_cupy, swir_dask_cupy, tir_dask_cupy)
-    general_output_checks(
-        red_dask_cupy, dask_cupy_result, ebbi_expected_results)
-
-
-def test_true_color_cpu():
+@pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
+def test_ebbi_gpu(red_data, swir1_data, tir_data, result_ebbi):
+    result = ebbi(red_data, swir1_data, tir_data)
+    general_output_checks(red_data, result, result_ebbi)
+
+
+# true_color ----------
+@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_true_color_numpy_equals_dask_numpy(random_data):
     # vanilla numpy version
-    red_numpy = create_test_arr(red_data)
-    green_numpy = create_test_arr(green_data)
-    blue_numpy = create_test_arr(blue_data)
+    red_numpy = create_test_raster(random_data, backend="numpy")
+    green_numpy = create_test_raster(random_data, backend="numpy")
+    blue_numpy = create_test_raster(random_data, backend="numpy")
     numpy_result = true_color(
         red_numpy, green_numpy, blue_numpy, name='np_true_color'
     )
@@ -582,14 +452,14 @@ def test_true_color_cpu():
     general_output_checks(red_numpy, numpy_result, verify_attrs=False)
 
     # dask
-    red_dask = create_test_arr(red_data, backend='dask')
-    green_dask = create_test_arr(green_data, backend='dask')
-    blue_dask = create_test_arr(blue_data, backend='dask')
+    red_dask = create_test_raster(random_data, backend='dask')
+    green_dask = create_test_raster(random_data, backend='dask')
+    blue_dask = create_test_raster(random_data, backend='dask')
     dask_result = true_color(
         red_dask, green_dask, blue_dask, name='dask_true_color'
     )
     assert dask_result.name == 'dask_true_color'
-    general_output_checks(red_numpy, numpy_result, verify_attrs=False)
+    general_output_checks(red_dask, dask_result, verify_attrs=False)
 
     np.testing.assert_allclose(
         numpy_result.data, dask_result.compute().data, equal_nan=True
diff --git a/xrspatial/tests/test_pathfinding.py b/xrspatial/tests/test_pathfinding.py
index 48e64fc..fd2de2e 100755
--- a/xrspatial/tests/test_pathfinding.py
+++ b/xrspatial/tests/test_pathfinding.py
@@ -1,31 +1,67 @@
 import numpy as np
-import xarray as xr
+import pytest
 
 from xrspatial import a_star_search
 
 from xrspatial.tests.general_checks import general_output_checks
+from xrspatial.tests.general_checks import create_test_raster
 
 
-def test_a_star_search():
-    agg = xr.DataArray(np.array([[0, 1, 0, 0],
-                                 [1, 1, 0, 0],
-                                 [0, 1, 2, 2],
-                                 [1, 0, 2, 0],
-                                 [0, 2, 2, 2]]),
-                       dims=['lat', 'lon'])
+@pytest.fixture
+def input_data():
+    data = np.array([[0, 1, 0, 0],
+                     [1, 1, 0, 0],
+                     [0, 1, 2, 2],
+                     [1, 0, 2, 0],
+                     [0, 2, 2, 2]])
+    agg = create_test_raster(data, dims=['lat', 'lon'])
+    return agg
 
-    height, width = agg.shape
-    _lon = np.linspace(0, width - 1, width)
-    _lat = np.linspace(height - 1, 0, height)
-    agg['lon'] = _lon
-    agg['lat'] = _lat
+
+@pytest.fixture
+def input_data_with_nans():
+    data = np.array([[0, 1, 0, 0],
+                     [1, 1, np.nan, 0],
+                     [0, 1, 2, 2],
+                     [1, 0, 2, 0],
+                     [0, np.nan, 2, 2]])
+    agg = create_test_raster(data, dims=['lat', 'lon'])
+
+    # start and end at a nan pixel, coordinate in (lat, lon) format
+    start = (1, 2)
+    goal = (4, 1)
+    return agg, start, goal
+
+
+@pytest.fixture
+def result_8_connectivity():
+    expected_result = np.array([[np.nan, np.nan, 0., np.nan],
+                                [np.nan, 1.41421356, np.nan, np.nan],
+                                [np.nan, 2.41421356, np.nan, np.nan],
+                                [np.nan, 3.41421356, np.nan, np.nan],
+                                [np.nan, np.nan, np.nan, np.nan]])
+    return expected_result
+
+
+@pytest.fixture
+def result_4_connectivity():
+    expected_result = np.array([[np.nan, 1, 0., np.nan],
+                                [np.nan, 2, np.nan, np.nan],
+                                [np.nan, 3, np.nan, np.nan],
+                                [np.nan, 4, np.nan, np.nan],
+                                [np.nan, np.nan, np.nan, np.nan]])
+    return expected_result
+
+
+def test_a_star_search_no_barriers(input_data):
+    agg = input_data
     barriers = []
     # no barriers, there always path from a start location to a goal location
-    for x0 in _lon:
-        for y0 in _lat:
+    for x0 in agg['lon']:
+        for y0 in agg['lat']:
             start = (y0, x0)
-            for x1 in _lon:
-                for y1 in _lat:
+            for x1 in agg['lon']:
+                for y1 in agg['lat']:
                     goal = (y1, x1)
                     path_agg = a_star_search(
                         agg, start, goal, barriers, 'lon', 'lat'
@@ -39,12 +75,15 @@ def test_a_star_search():
                         assert np.nanmax(path_agg) > 0
                         assert np.nanmin(path_agg) == 0
 
+
+def test_a_star_search_with_barriers(input_data):
+    agg = input_data
     barriers = [1]
     # set pixels with value 1 as barriers,
     # cannot go from (0, 0) to anywhere since it is surrounded by 1s
-    start = (4, 0)
-    for x1 in _lon:
-        for y1 in _lat:
+    start = (0, 0)
+    for x1 in agg['lon']:
+        for y1 in agg['lat']:
             goal = (y1, x1)
             if (goal != start):
                 path_agg = a_star_search(
@@ -54,62 +93,46 @@ def test_a_star_search():
                 expected_results = np.full(agg.shape, np.nan)
                 general_output_checks(agg, path_agg, expected_results)
 
-    # test with nans
-    agg = xr.DataArray(np.array([[0, 1, 0, 0],
-                                 [1, 1, np.nan, 0],
-                                 [0, 1, 2, 2],
-                                 [1, 0, 2, 0],
-                                 [0, np.nan, 2, 2]]),
-                       dims=['lat', 'lon'])
-
-    height, width = agg.shape
-    _lon = np.linspace(0, width - 1, width)
-    _lat = np.linspace(0, height - 1, height)
-    agg['lon'] = _lon
-    agg['lat'] = _lat
-    # start and end at a nan pixel, coordinate in (lat, lon) format
-    start = (1, 2)
-    goal = (4, 1)
+
+def test_a_star_search_snap(input_data_with_nans):
+    agg, start, goal = input_data_with_nans
+
     # no barriers
     barriers = []
     # no snap
     no_snap_path_agg = a_star_search(agg, start, goal, barriers, 'lon', 'lat')
     # no path, all cells in path_agg are nans
-    assert np.isnan(no_snap_path_agg).all()
+    np.testing.assert_array_equal(no_snap_path_agg, np.nan)
 
     # set snap_start = True, snap_goal = False
-    snap_start_path_agg = a_star_search(agg, start, goal, barriers,
-                                        'lon', 'lat', snap_start=True)
+    snap_start_path_agg = a_star_search(agg, start, goal, barriers, 'lon', 'lat', snap_start=True)
     # no path, all cells in path_agg are nans
-    assert np.isnan(snap_start_path_agg).all()
+    np.testing.assert_array_equal(snap_start_path_agg, np.nan)
 
     # set snap_start = False, snap_goal = True
-    snap_goal_path_agg = a_star_search(agg, start, goal, barriers,
-                                       'lon', 'lat', snap_goal=True)
+    snap_goal_path_agg = a_star_search(agg, start, goal, barriers, 'lon', 'lat', snap_goal=True)
     # no path, all cells in path_agg are nans
-    assert np.isnan(snap_goal_path_agg).all()
+    np.testing.assert_array_equal(snap_goal_path_agg, np.nan)
+
+
+def test_a_star_search_connectivity(
+    input_data_with_nans,
+    result_8_connectivity,
+    result_4_connectivity
+):
+    agg, start, goal = input_data_with_nans
+    # no barriers
+    barriers = []
 
     # set snap_start = True, snap_goal = True
     # 8-connectivity as default
-    path_agg_8 = a_star_search(agg, start, goal, barriers, 'lon', 'lat',
-                               snap_start=True, snap_goal=True)
-    # path exists
-    expected_result_8 = np.array([[np.nan, np.nan, 0., np.nan],
-                                 [np.nan, 1.41421356, np.nan, np.nan],
-                                 [np.nan, 2.41421356, np.nan, np.nan],
-                                 [np.nan, 3.41421356, np.nan, np.nan],
-                                 [np.nan, np.nan, np.nan, np.nan]])
-    assert ((np.isnan(path_agg_8) & np.isnan(expected_result_8)) | (
-                    abs(path_agg_8 - expected_result_8) <= 1e-5)).all()
+    path_agg_8 = a_star_search(
+        agg, start, goal, barriers, 'lon', 'lat', snap_start=True, snap_goal=True
+    )
+    np.testing.assert_allclose(path_agg_8, result_8_connectivity, equal_nan=True)
 
     # 4-connectivity
-    path_agg_4 = a_star_search(agg, start, goal, barriers, 'lon', 'lat',
-                               snap_start=True, snap_goal=True, connectivity=4)
-    # path exists
-    expected_result_4 = np.array([[np.nan, 1, 0., np.nan],
-                                  [np.nan, 2, np.nan, np.nan],
-                                  [np.nan, 3, np.nan, np.nan],
-                                  [np.nan, 4, np.nan, np.nan],
-                                  [np.nan, np.nan, np.nan, np.nan]])
-    assert ((np.isnan(path_agg_4) & np.isnan(expected_result_4)) | (
-                    abs(path_agg_4 - expected_result_4) <= 1e-5)).all()
+    path_agg_4 = a_star_search(
+        agg, start, goal, barriers, 'lon', 'lat', snap_start=True, snap_goal=True, connectivity=4
+    )
+    np.testing.assert_allclose(path_agg_4, result_4_connectivity, equal_nan=True)
diff --git a/xrspatial/tests/test_proximity.py b/xrspatial/tests/test_proximity.py
index 5c5b57e..c6c0cfc 100755
--- a/xrspatial/tests/test_proximity.py
+++ b/xrspatial/tests/test_proximity.py
@@ -22,172 +22,198 @@ def test_great_circle_distance():
             assert e_info
 
 
-def create_test_raster():
+@pytest.fixture
+def test_raster(backend):
     height, width = 4, 6
+    # create test raster, all non-zero cells are unique,
+    # this is to test allocation and direction against corresponding proximity
     data = np.asarray([[0., 0., 0., 0., 0., 2.],
                        [0., 0., 1., 0., 0., 0.],
                        [0., np.inf, 3., 0., 0., 0.],
                        [4., 0., 0., 0., np.nan, 0.]])
     _lon = np.linspace(-20, 20, width)
     _lat = np.linspace(20, -20, height)
-
-    numpy_agg = xr.DataArray(data, dims=['lat', 'lon'])
-    dask_numpy_agg = xr.DataArray(
-        da.from_array(data, chunks=(4, 3)), dims=['lat', 'lon'])
-    numpy_agg['lon'] = dask_numpy_agg['lon'] = _lon
-    numpy_agg['lat'] = dask_numpy_agg['lat'] = _lat
-
-    return numpy_agg, dask_numpy_agg
+    raster = xr.DataArray(data, dims=['lat', 'lon'])
+    raster['lon'] = _lon
+    raster['lat'] = _lat
+    if 'dask' in backend:
+        raster.data = da.from_array(data, chunks=(4, 3))
+    return raster
 
 
-def test_proximity():
-    raster_numpy, raster_dask = create_test_raster()
-
+@pytest.fixture
+def result_default_proximity():
     # DEFAULT SETTINGS
-    expected_results_default = np.array([
+    expected_result = np.array([
         [20.82733247, 15.54920505, 13.33333333, 15.54920505,  8., 0.],
         [16., 8., 0., 8., 15.54920505, 13.33333333],
         [13.33333333, 8., 0., 8., 16., 24.],
         [0., 8., 13.33333333, 15.54920505, 20.82733247, 27.45501371]
-    ])
-    # numpy case
-    default_prox = proximity(raster_numpy, x='lon', y='lat')
-    general_output_checks(raster_numpy, default_prox, expected_results_default)
-    # dask case
-    default_prox_dask = proximity(raster_dask, x='lon', y='lat')
-    general_output_checks(
-        raster_dask, default_prox_dask, expected_results_default)
+    ], dtype=np.float32)
+    return expected_result
 
-    # TARGET VALUES SETTING
+
+@pytest.fixture
+def result_target_proximity():
     target_values = [2, 3]
-    expected_results_target = np.array([
+    expected_result = np.array([
         [31.09841011, 27.84081736, 24., 16., 8., 0.],
-        [20.82733247, 15.54920505, 13.33333333, 15.54920505, 15.54920505, 13.33333333],  # noqa
+        [20.82733247, 15.54920505, 13.33333333, 15.54920505, 15.54920505, 13.33333333],
         [16., 8., 0., 8., 16., 24.],
-        [20.82733247, 15.54920505, 13.33333333, 15.54920505, 20.82733247, 27.45501371]  # noqa
-    ])
-    # numpy case
-    target_prox = proximity(
-        raster_numpy, x='lon', y='lat', target_values=target_values)
-    general_output_checks(raster_numpy, target_prox, expected_results_target)
-    # dask case
-    target_prox_dask = proximity(
-        raster_dask, x='lon', y='lat', target_values=target_values)
-    general_output_checks(
-        raster_dask, target_prox_dask, expected_results_target)
+        [20.82733247, 15.54920505, 13.33333333, 15.54920505, 20.82733247, 27.45501371]
+    ], dtype=np.float32)
+    return target_values, expected_result
 
+
+@pytest.fixture
+def result_manhattan_proximity():
     # distance_metric SETTING: MANHATTAN
-    expected_results_manhattan = np.array([
+    expected_result = np.array([
         [29.33333333, 21.33333333, 13.33333333, 16., 8., 0.],
         [16., 8., 0., 8., 16., 13.33333333],
         [13.33333333, 8., 0., 8., 16., 24.],
         [0., 8., 13.33333333, 21.33333333, 29.33333333, 37.33333333]
-    ])
-    # numpy case
-    manhattan_prox = proximity(
-        raster_numpy, x='lon', y='lat', distance_metric='MANHATTAN')
-    general_output_checks(
-        raster_numpy, manhattan_prox, expected_results_manhattan)
-    # dask case
-    manhattan_prox_dask = proximity(raster_dask, x='lon', y='lat',
-                                    distance_metric='MANHATTAN')
-    general_output_checks(
-        raster_dask, manhattan_prox_dask, expected_results_manhattan)
+    ], dtype=np.float32)
+    return expected_result
+
 
+@pytest.fixture
+def result_great_circle_proximity():
     # distance_metric SETTING: GREAT_CIRCLE
-    expected_results_great_circle = np.array([
-        [2278099.27025501, 1717528.97437217, 1484259.87724365, 1673057.17235307, 836769.1780019, 0.],  # noqa
-        [1768990.54084204, 884524.60324856, 0., 884524.60324856, 1717528.97437217, 1484259.87724365],  # noqa
-        [1484259.87724365, 884524.60324856, 0., 884524.60324856, 1768990.54084204, 2653336.85436932],  # noqa
-        [0., 836769.1780019, 1484259.87724365, 1717528.97437217, 2278099.27025501, 2986647.12982316]  # noqa
-    ])
-    great_circle_prox = proximity(raster_numpy, x='lon', y='lat',
-                                  distance_metric='GREAT_CIRCLE')
-    general_output_checks(
-        raster_numpy, great_circle_prox, expected_results_great_circle)
-    # dask case
-    great_circle_prox_dask = proximity(
-        raster_dask, x='lon', y='lat', distance_metric='GREAT_CIRCLE'
-    )
-    general_output_checks(
-        raster_dask, great_circle_prox_dask, expected_results_great_circle)
+    expected_result = np.array([
+        [2278099.27025501, 1717528.97437217, 1484259.87724365, 1673057.17235307, 836769.1780019, 0],
+        [1768990.54084204, 884524.60324856, 0, 884524.60324856, 1717528.97437217, 1484259.87724365],
+        [1484259.87724365, 884524.60324856, 0, 884524.60324856, 1768990.54084204, 2653336.85436932],
+        [0, 836769.1780019, 1484259.87724365, 1717528.97437217, 2278099.27025501, 2986647.12982316]
+    ], dtype=np.float32)
+    return expected_result
 
+
+@pytest.fixture
+def result_max_distance_proximity():
     # max_distance setting
     max_distance = 10
-    expected_result_max_distance = np.array([
+    expected_result = np.array([
         [np.nan, np.nan, np.nan, np.nan, 8., 0.],
         [np.nan, 8., 0., 8., np.nan, np.nan],
         [np.nan, 8., 0., 8., np.nan, np.nan],
         [0., 8., np.nan, np.nan, np.nan, np.nan]
-    ])
-    # numpy case
-    max_distance_prox = proximity(
-        raster_numpy, x='lon', y='lat', max_distance=max_distance
-    )
-    general_output_checks(
-        raster_numpy, max_distance_prox, expected_result_max_distance)
-    # dask case
-    max_distance_prox_dask = proximity(
-        raster_dask, x='lon', y='lat', max_distance=max_distance
-    )
-    general_output_checks(
-        raster_dask, max_distance_prox_dask, expected_result_max_distance
-    )
+    ], dtype=np.float32)
+    return max_distance, expected_result
 
 
-def test_allocation():
-    # create test raster, all non-zero cells are unique,
-    # this is to test against corresponding proximity
-    raster_numpy, raster_dask = create_test_raster()
-    expected_results = np.array([
+@pytest.fixture
+def result_default_allocation():
+    expected_result = np.array([
         [1., 1., 1., 1., 2., 2.],
         [1., 1., 1., 1., 2., 2.],
         [4., 3., 3., 3., 3., 3.],
         [4., 4., 3., 3., 3., 3.]
-    ])
-    allocation_agg = allocation(raster_numpy, x='lon', y='lat')
-    general_output_checks(raster_numpy, allocation_agg, expected_results)
-
-    # check against corresponding proximity
-    proximity_agg = proximity(raster_numpy, x='lon', y='lat')
-    xcoords = allocation_agg['lon'].data
-    ycoords = allocation_agg['lat'].data
-    for y in range(raster_numpy.shape[0]):
-        for x in range(raster_numpy.shape[1]):
-            a = allocation_agg.data[y, x]
-            py, px = np.where(raster_numpy.data == a)
-            # non-zero cells in raster are unique, thus len(px)=len(py)=1
-            d = euclidean_distance(xcoords[x], xcoords[px[0]],
-                                   ycoords[y], ycoords[py[0]])
-            assert proximity_agg.data[y, x] == np.float32(d)
+    ], dtype=np.float32)
+    return expected_result
 
-    # dask case
-    allocation_agg_dask = allocation(raster_dask, x='lon', y='lat')
-    general_output_checks(raster_dask, allocation_agg_dask, expected_results)
 
+@pytest.fixture
+def result_max_distance_allocation():
     # max_distance setting
     max_distance = 10
-    expected_results_max_distance = np.array([
+    expected_result = np.array([
         [np.nan, np.nan, np.nan, np.nan, 2., 2.],
         [np.nan, 1., 1., 1., np.nan, np.nan],
         [np.nan, 3., 3., 3., np.nan, np.nan],
         [4., 4., np.nan, np.nan, np.nan, np.nan]
-    ])
-    # numpy case
-    max_distance_alloc = allocation(
-        raster_numpy, x='lon', y='lat', max_distance=max_distance
-    )
+    ], dtype=np.float32)
+    return max_distance, expected_result
+
+
+@pytest.fixture
+def result_default_direction():
+    expected_result = np.array([
+        [50.194427, 30.963757, 360., 329.03625, 90., 0.],
+        [90., 90., 0., 270., 149.03624, 180.],
+        [360., 90., 0., 270., 270., 270.],
+        [0., 270., 180., 210.96376, 230.19443, 240.9454]
+    ], dtype=np.float32)
+    return expected_result
+
+
+@pytest.fixture
+def result_max_distance_direction():
+    # max_distance setting
+    max_distance = 10
+    expected_result = np.array([
+        [np.nan, np.nan, np.nan, np.nan, 90., 0.],
+        [np.nan, 90., 0., 270., np.nan, np.nan],
+        [np.nan, 90., 0., 270., np.nan, np.nan],
+        [0., 270., np.nan, np.nan, np.nan, np.nan]
+    ], dtype=np.float32)
+    return max_distance, expected_result
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_default_proximity(test_raster, result_default_proximity):
+    default_prox = proximity(test_raster, x='lon', y='lat')
+    general_output_checks(test_raster, default_prox, result_default_proximity, verify_dtype=True)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_target_proximity(test_raster, result_target_proximity):
+    target_values, expected_result = result_target_proximity
+    target_prox = proximity(test_raster, x='lon', y='lat', target_values=target_values)
+    general_output_checks(test_raster, target_prox, expected_result, verify_dtype=True)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_manhattan_proximity(test_raster, result_manhattan_proximity):
+    manhattan_prox = proximity(test_raster, x='lon', y='lat', distance_metric='MANHATTAN')
     general_output_checks(
-        raster_numpy, max_distance_alloc, expected_results_max_distance)
-    # dask case
-    max_distance_alloc_dask = allocation(
-        raster_dask, x='lon', y='lat', max_distance=max_distance
+        test_raster, manhattan_prox, result_manhattan_proximity, verify_dtype=True
     )
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_great_circle_proximity(test_raster, result_great_circle_proximity):
+    great_circle_prox = proximity(test_raster, x='lon', y='lat', distance_metric='GREAT_CIRCLE')
     general_output_checks(
-        raster_dask, max_distance_alloc_dask, expected_results_max_distance
+        test_raster, great_circle_prox, result_great_circle_proximity, verify_dtype=True
     )
 
 
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_max_distance_proximity(test_raster, result_max_distance_proximity):
+    max_distance, expected_result = result_max_distance_proximity
+    max_distance_prox = proximity(test_raster, x='lon', y='lat', max_distance=max_distance)
+    general_output_checks(test_raster, max_distance_prox, expected_result, verify_dtype=True)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_default_allocation(test_raster, result_default_allocation):
+    allocation_agg = allocation(test_raster, x='lon', y='lat')
+    general_output_checks(test_raster, allocation_agg, result_default_allocation, verify_dtype=True)
+
+
+@pytest.mark.parametrize("backend", ['numpy'])
+def test_default_allocation_against_proximity(test_raster, result_default_proximity):
+    allocation_agg = allocation(test_raster, x='lon', y='lat')
+    # check against corresponding proximity
+    xcoords = allocation_agg['lon'].data
+    ycoords = allocation_agg['lat'].data
+    for y in range(test_raster.shape[0]):
+        for x in range(test_raster.shape[1]):
+            a = allocation_agg.data[y, x]
+            py, px = np.where(test_raster.data == a)
+            # non-zero cells in raster are unique, thus len(px)=len(py)=1
+            d = euclidean_distance(xcoords[x], xcoords[px[0]], ycoords[y], ycoords[py[0]])
+            np.testing.assert_allclose(result_default_proximity[y, x], d)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_max_distance_allocation(test_raster, result_max_distance_allocation):
+    max_distance, expected_result = result_max_distance_allocation
+    max_distance_alloc = allocation(test_raster, x='lon', y='lat', max_distance=max_distance)
+    general_output_checks(test_raster, max_distance_alloc, expected_result, verify_dtype=True)
+
+
 def test_calc_direction():
     n = 3
     x1, y1 = 1, 1
@@ -204,57 +230,28 @@ def test_calc_direction():
     assert (abs(output-expected_output) <= tolerance).all()
 
 
-def test_direction():
-    raster_numpy, raster_dask = create_test_raster()
-
-    expected_results = np.array([
-        [50.194427, 30.963757, 360., 329.03625, 90., 0.],
-        [90., 90., 0., 270., 149.03624, 180.],
-        [360., 90., 0., 270., 270., 270.],
-        [0., 270., 180., 210.96376, 230.19443, 240.9454]
-    ], dtype=np.float32)
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_default_direction(test_raster, result_default_direction):
+    direction_agg = direction(test_raster, x='lon', y='lat')
+    general_output_checks(test_raster, direction_agg, result_default_direction)
 
-    # numpy case
-    direction_agg = direction(raster_numpy, x='lon', y='lat')
-    # output must be an xarray DataArray
-    general_output_checks(raster_numpy, direction_agg, expected_results)
 
-    # test against allocation
-    allocation_agg = allocation(raster_numpy, x='lon', y='lat')
-    xcoords = allocation_agg['lon'].data
-    ycoords = allocation_agg['lat'].data
-    for y in range(raster_numpy.shape[0]):
-        for x in range(raster_numpy.shape[1]):
-            a = allocation_agg.data[y, x]
-            py, px = np.where(raster_numpy.data == a)
+@pytest.mark.parametrize("backend", ['numpy'])
+def test_default_direction_against_allocation(test_raster, result_default_allocation):
+    direction_agg = direction(test_raster, x='lon', y='lat')
+    xcoords = direction_agg['lon'].data
+    ycoords = direction_agg['lat'].data
+    for y in range(test_raster.shape[0]):
+        for x in range(test_raster.shape[1]):
+            a = result_default_allocation.data[y, x]
+            py, px = np.where(test_raster.data == a)
             # non-zero cells in raster are unique, thus len(px)=len(py)=1
-            d = _calc_direction(xcoords[x], xcoords[px[0]],
-                                ycoords[y], ycoords[py[0]])
-            assert direction_agg.data[y, x] == d
+            d = _calc_direction(xcoords[x], xcoords[px[0]], ycoords[y], ycoords[py[0]])
+            np.testing.assert_allclose(direction_agg.data[y, x], d)
 
-    # dask case
-    direction_agg_dask = direction(raster_dask, x='lon', y='lat')
-    general_output_checks(raster_dask, direction_agg_dask, expected_results)
 
-    # max_distance setting
-    max_distance = 10
-    expected_results_max_distance = np.array([
-        [np.nan, np.nan, np.nan, np.nan, 90., 0.],
-        [np.nan, 90., 0., 270., np.nan, np.nan],
-        [np.nan, 90., 0., 270., np.nan, np.nan],
-        [0., 270., np.nan, np.nan, np.nan, np.nan]
-    ], dtype=np.float32)
-    # numpy case
-    max_distance_direction = direction(
-        raster_numpy, x='lon', y='lat', max_distance=max_distance
-    )
-    general_output_checks(
-        raster_numpy, max_distance_direction, expected_results_max_distance
-    )
-    # dask case
-    max_distance_direction_dask = direction(
-        raster_dask, x='lon', y='lat', max_distance=max_distance
-    )
-    general_output_checks(
-        raster_dask, max_distance_direction_dask, expected_results_max_distance
-    )
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_max_distance_direction(test_raster, result_max_distance_direction):
+    max_distance, expected_result = result_max_distance_direction
+    max_distance_direction = direction(test_raster, x='lon', y='lat', max_distance=max_distance)
+    general_output_checks(test_raster, max_distance_direction, expected_result, verify_dtype=True)
diff --git a/xrspatial/tests/test_slope.py b/xrspatial/tests/test_slope.py
index cd9b1a5..225205d 100755
--- a/xrspatial/tests/test_slope.py
+++ b/xrspatial/tests/test_slope.py
@@ -1,112 +1,76 @@
 import pytest
-import xarray as xr
 import numpy as np
 
-import dask.array as da
-
 from xrspatial import slope
 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
 
 
-# Test Data -----------------------------------------------------------------
-
-'''
-Notes:
-------
-The `elevation` data was run through QGIS slope function to
-get values to compare against.  Xarray-Spatial currently handles
-edges by padding with nan which is different than QGIS but acknowledged
-'''
-
-elevation = np.asarray([
-    [1432.6542, 1432.4764, 1432.4764, 1432.1207, 1431.9429, np.nan],
-    [1432.6542, 1432.6542, 1432.4764, 1432.2986, 1432.1207, np.nan],
-    [1432.832, 1432.6542, 1432.4764, 1432.2986, 1432.1207, np.nan],
-    [1432.832, 1432.6542, 1432.4764, 1432.4764, 1432.1207, np.nan],
-    [1432.832, 1432.6542, 1432.6542, 1432.4764, 1432.2986, np.nan],
-    [1432.832, 1432.6542, 1432.6542, 1432.4764, 1432.2986, np.nan],
-    [1432.832, 1432.832, 1432.6542, 1432.4764, 1432.4764, np.nan]],
-    dtype=np.float32
-)
+def input_data(backend):
+    # Notes:
+    # ------
+    # The `elevation` data was run through QGIS slope function to
+    # get values to compare against.  Xarray-Spatial currently handles
+    # edges by padding with nan which is different than QGIS but acknowledged
+
+    elevation = np.asarray([
+        [1432.6542, 1432.4764, 1432.4764, 1432.1207, 1431.9429, np.nan],
+        [1432.6542, 1432.6542, 1432.4764, 1432.2986, 1432.1207, np.nan],
+        [1432.832, 1432.6542, 1432.4764, 1432.2986, 1432.1207, np.nan],
+        [1432.832, 1432.6542, 1432.4764, 1432.4764, 1432.1207, np.nan],
+        [1432.832, 1432.6542, 1432.6542, 1432.4764, 1432.2986, np.nan],
+        [1432.832, 1432.6542, 1432.6542, 1432.4764, 1432.2986, np.nan],
+        [1432.832, 1432.832, 1432.6542, 1432.4764, 1432.4764, np.nan]],
+        dtype=np.float32
+    )
+    raster = create_test_raster(elevation, backend, attrs={'res': (10, 10)})
+    return raster
 
-qgis_slope = np.asarray(
-    [[0.8052942, 0.742317, 1.1390567, 1.3716657, np.nan, np.nan],
-     [0.74258685, 0.742317, 1.0500116, 1.2082565, np.nan, np.nan],
-     [0.56964326, 0.9002944, 0.9002944, 1.0502871, np.nan, np.nan],
-     [0.5095078, 0.9003686, 0.742317, 1.1390567, np.nan, np.nan],
-     [0.6494868, 0.64938396, 0.5692523, 1.0500116, np.nan, np.nan],
-     [0.80557066, 0.56964326, 0.64914393, 0.9002944, np.nan, np.nan],
-     [0.6494868, 0.56964326, 0.8052942, 0.742317, np.nan, np.nan]],
-    dtype=np.float32)
 
+@pytest.fixture
+def qgis_output():
+    qgis_slope = np.asarray(
+        [[0.8052942, 0.742317, 1.1390567, 1.3716657, np.nan, np.nan],
+         [0.74258685, 0.742317, 1.0500116, 1.2082565, np.nan, np.nan],
+         [0.56964326, 0.9002944, 0.9002944, 1.0502871, np.nan, np.nan],
+         [0.5095078, 0.9003686, 0.742317, 1.1390567, np.nan, np.nan],
+         [0.6494868, 0.64938396, 0.5692523, 1.0500116, np.nan, np.nan],
+         [0.80557066, 0.56964326, 0.64914393, 0.9002944, np.nan, np.nan],
+         [0.6494868, 0.56964326, 0.8052942, 0.742317, np.nan, np.nan]],
+        dtype=np.float32)
+    return qgis_slope
 
-def test_slope_against_qgis_cpu():
 
+def test_numpy_equals_qgis(qgis_output):
     # slope by xrspatial
-    agg_numpy = xr.DataArray(elevation, attrs={'res': (10.0, 10.0)})
-    xrspatial_slope_numpy = slope(agg_numpy, name='slope_numpy')
-    general_output_checks(agg_numpy, xrspatial_slope_numpy)
+    numpy_agg = input_data(backend='numpy')
+    xrspatial_slope_numpy = slope(numpy_agg, name='slope_numpy')
+    general_output_checks(numpy_agg, xrspatial_slope_numpy)
     assert xrspatial_slope_numpy.name == 'slope_numpy'
 
-    agg_dask = xr.DataArray(
-        da.from_array(elevation, chunks=(3, 3)), attrs={'res': (10.0, 10.0)})
-    xrspatial_slope_dask = slope(agg_dask, name='slope_dask')
-    general_output_checks(agg_dask, xrspatial_slope_dask)
-    assert xrspatial_slope_dask.name == 'slope_dask'
-
-    # numpy and dask case produce same results
-    np.testing.assert_allclose(
-        xrspatial_slope_numpy.data,
-        xrspatial_slope_dask.compute().data,
-        equal_nan=True
-    )
-
-    # nan border edges
-    xrspatial_edges = [
-        xrspatial_slope_numpy.data[0, :],
-        xrspatial_slope_numpy.data[-1, :],
-        xrspatial_slope_numpy.data[:, 0],
-        xrspatial_slope_numpy.data[:, -1],
-    ]
-    for edge in xrspatial_edges:
-        np.testing.assert_allclose(
-            edge, np.full(edge.shape, np.nan), equal_nan=True
-        )
-
-    # test against QGIS
-    xrspatial_vals = xrspatial_slope_numpy.values[1:-1, 1:-1]
-    qgis_vals = qgis_slope[1:-1, 1:-1]
+    xrspatial_vals = xrspatial_slope_numpy.data[1:-1, 1:-1]
+    qgis_vals = qgis_output[1:-1, 1:-1]
     np.testing.assert_allclose(xrspatial_vals, qgis_vals, equal_nan=True)
 
+    # nan border edges
+    assert_nan_edges_effect(xrspatial_slope_numpy)
 
-@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
-def test_slope_against_qgis_gpu():
 
-    import cupy
+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, slope)
 
-    # slope by xrspatial
-    agg_cupy = xr.DataArray(
-        cupy.asarray(elevation), attrs={'res': (10.0, 10.0)})
-    xrspatial_slope_gpu = slope(agg_cupy, name='slope_cupy')
-    general_output_checks(agg_cupy, xrspatial_slope_gpu)
-    assert xrspatial_slope_gpu.name == 'slope_cupy'
 
-    agg_numpy = xr.DataArray(elevation, attrs={'res': (10.0, 10.0)})
-    xrspatial_slope_cpu = slope(agg_numpy, name='slope_numpy')
-
-    # both cpu and gpu produce same results
-    np.testing.assert_allclose(
-        xrspatial_slope_cpu.data,
-        xrspatial_slope_gpu.data.get(),
-        equal_nan=True
-    )
-
-    # test against QGIS
-    # nan border edges
-    xrspatial_vals = xrspatial_slope_gpu.data[1:-1, 1:-1].get()
-    qgis_vals = qgis_slope[1:-1, 1:-1]
-    np.testing.assert_allclose(
-        xrspatial_vals, qgis_vals, equal_nan=True
-    )
+@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('numpy')
+    cupy_agg = input_data('cupy')
+    assert_numpy_equals_cupy(numpy_agg, cupy_agg, slope)
diff --git a/xrspatial/tests/test_terrain.py b/xrspatial/tests/test_terrain.py
index 476d561..9651a99 100755
--- a/xrspatial/tests/test_terrain.py
+++ b/xrspatial/tests/test_terrain.py
@@ -36,10 +36,7 @@ def test_terrain_cpu():
     assert isinstance(terrain_dask.data, da.Array)
 
     terrain_dask = terrain_dask.compute()
-    assert np.isclose(
-        terrain_numpy.data, terrain_dask.data,
-        rtol=1e-05, atol=1e-07, equal_nan=True
-    ).all()
+    np.testing.assert_allclose(terrain_numpy.data, terrain_dask.data, rtol=1e-05, atol=1e-07)
 
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
@@ -52,7 +49,4 @@ def test_terrain_gpu():
     data_cupy = create_test_arr(backend='cupy')
     terrain_cupy = generate_terrain(data_cupy)
 
-    assert np.isclose(
-        terrain_numpy.data, terrain_cupy.data,
-        rtol=1e-05, atol=1e-07, equal_nan=True
-    ).all()
+    np.testing.assert_allclose(terrain_numpy.data, terrain_cupy.data.get(), rtol=1e-05, atol=1e-07)
diff --git a/xrspatial/tests/test_viewshed.py b/xrspatial/tests/test_viewshed.py
index 9247d35..459f7d6 100644
--- a/xrspatial/tests/test_viewshed.py
+++ b/xrspatial/tests/test_viewshed.py
@@ -1,67 +1,60 @@
-import xarray as xa
 import pytest
 
-import datashader as ds
-from ..gpu_rtx import has_rtx
-from xrspatial import viewshed
-
 import numpy as np
 import pandas as pd
+import xarray as xa
 
-# create an empty image of size 5*5
-H = 5
-W = 5
+import datashader as ds
+from ..gpu_rtx import has_rtx
+from xrspatial import viewshed
 
-canvas = ds.Canvas(plot_width=W, plot_height=H,
-                   x_range=(-20, 20), y_range=(-20, 20))
+from xrspatial.tests.general_checks import general_output_checks
 
-empty_df = pd.DataFrame({
-   'x': np.random.normal(.5, 1, 0),
-   'y': np.random.normal(.5, 1, 0)
-})
-empty_agg = canvas.points(empty_df, 'x', 'y')
 
-# coordinates
-xs = empty_agg.coords['x'].values
-ys = empty_agg.coords['x'].values
+@pytest.fixture
+def empty_agg():
+    # create an empty image of size 5*5
+    H = 5
+    W = 5
 
-# define some values for observer's elevation to test
-OBS_ELEVS = [-1, 0, 1]
+    canvas = ds.Canvas(plot_width=W, plot_height=H,
+                       x_range=(-20, 20), y_range=(-20, 20))
 
-TERRAIN_ELEV_AT_VP = [-1, 0, 1]
+    empty_df = pd.DataFrame({
+       'x': np.random.normal(.5, 1, 0),
+       'y': np.random.normal(.5, 1, 0)
+    })
+    agg = canvas.points(empty_df, 'x', 'y')
+    return agg
 
 
-def test_viewshed_invalid_x_view():
+def test_viewshed_invalid_x_view(empty_agg):
+    xs = empty_agg.coords['x'].values
     OBSERVER_X = xs[0] - 1
     OBSERVER_Y = 0
-    with pytest.raises(Exception) as e_info:  # NOQA
-        viewshed(raster=empty_agg, x=OBSERVER_X, y=OBSERVER_Y,
-                 observer_elev=10)
+    with pytest.raises(Exception):
+        viewshed(raster=empty_agg, x=OBSERVER_X, y=OBSERVER_Y, observer_elev=10)
 
 
-def test_viewshed_invalid_y_view():
+def test_viewshed_invalid_y_view(empty_agg):
+    ys = empty_agg.coords['y'].values
     OBSERVER_X = 0
     OBSERVER_Y = ys[-1] + 1
-    with pytest.raises(Exception) as e_info:  # NOQA
-        viewshed(raster=empty_agg, x=OBSERVER_X, y=OBSERVER_Y,
-                 observer_elev=10)
+    with pytest.raises(Exception):
+        viewshed(raster=empty_agg, x=OBSERVER_X, y=OBSERVER_Y, observer_elev=10)
 
 
-def test_viewshed_output_properties():
-    for obs_elev in OBS_ELEVS:
-        OBSERVER_X = xs[0]
-        OBSERVER_Y = ys[0]
-        v = viewshed(raster=empty_agg, x=OBSERVER_X, y=OBSERVER_Y,
-                     observer_elev=obs_elev)
+def test_viewshed(empty_agg):
+    H, W = empty_agg.shape
 
-        assert v.shape[0] == empty_agg.shape[0]
-        assert v.shape[1] == empty_agg.shape[1]
-        assert isinstance(v, xa.DataArray)
-        assert isinstance(v.values, np.ndarray)
-        assert type(v.values[0, 0]) == np.float64
+    # coordinates
+    xs = empty_agg.coords['x'].values
+    ys = empty_agg.coords['y'].values
 
+    # define some values for observer's elevation to test
+    OBS_ELEVS = [-1, 0, 1]
+    TERRAIN_ELEV_AT_VP = [-1, 0, 1]
 
-def test_viewshed():
     # check if a matrix is symmetric
     def check_symmetric(matrix, rtol=1e-05, atol=1e-08):
         return np.allclose(matrix, matrix.T, rtol=rtol, atol=atol)
@@ -85,8 +78,8 @@ def test_viewshed():
         tmp = min(height - y, x)
         s_bottom_y, s_left_x = y + tmp, x - tmp
 
-        return (f_top_y, f_left_x, f_bottom_y + 1, f_right_x + 1), \
-               (s_top_y, s_left_x, s_bottom_y + 1, s_right_x + 1)
+        return ((f_top_y, f_left_x, f_bottom_y + 1, f_right_x + 1),
+                (s_top_y, s_left_x, s_bottom_y + 1, s_right_x + 1))
 
     # test on 3 scenarios:
     #   empty image.
@@ -107,14 +100,17 @@ def test_viewshed():
         for elev_at_vp in TERRAIN_ELEV_AT_VP:
             for col_id, x in enumerate(xs):
                 for row_id, y in enumerate(ys):
+
                     empty_agg.values[row_id, col_id] = elev_at_vp
-                    v = viewshed(raster=empty_agg, x=x, y=y,
-                                 observer_elev=obs_elev)
+                    v = viewshed(raster=empty_agg, x=x, y=y, observer_elev=obs_elev)
+
+                    # validate output properties
+                    general_output_checks(empty_agg, v)
+
                     # angle at viewpoint is always 180
                     assert v[row_id, col_id] == 180
 
-                    if obs_elev + elev_at_vp >= 0 and \
-                            obs_elev >= abs(elev_at_vp):
+                    if obs_elev + elev_at_vp >= 0 and obs_elev >= abs(elev_at_vp):
                         # all cells are visible
                         assert (v.values > -1).all()
 
diff --git a/xrspatial/tests/test_zonal.py b/xrspatial/tests/test_zonal.py
index 1c2511d..ed3436d 100755
--- a/xrspatial/tests/test_zonal.py
+++ b/xrspatial/tests/test_zonal.py
@@ -1,3 +1,5 @@
+import pytest
+
 import numpy as np
 import pandas as pd
 import xarray as xr
@@ -10,65 +12,45 @@ from xrspatial import zonal_crosstab as crosstab
 from xrspatial import suggest_zonal_canvas
 from xrspatial import trim
 from xrspatial import crop
+from xrspatial.zonal import regions
 
+from xrspatial.tests.general_checks import create_test_raster
 
-from xrspatial.zonal import regions
 
+@pytest.fixture
+def data_zones(backend):
+    data = np.array([[0, 0, 1, 1, 2, 2, 3, 3],
+                     [0, 0, 1, 1, 2, 2, 3, 3],
+                     [0, 0, 1, 1, 2, np.nan, 3, 3]])
+    agg = create_test_raster(data, backend)
+    return agg
 
-def create_zones_values(backend):
-    zones_val = np.array([[0, 0, 1, 1, 2, 2, 3, 3],
-                          [0, 0, 1, 1, 2, 2, 3, 3],
-                          [0, 0, 1, 1, 2, np.nan, 3, 3]])
-    zones = xr.DataArray(zones_val)
 
-    values_val_2d = np.asarray([
+@pytest.fixture
+def data_values_2d(backend):
+    data = np.asarray([
         [0, 0, 1, 1, 2, 2, 3, np.inf],
         [0, 0, 1, 1, 2, np.nan, 3, 0],
         [np.inf, 0, 1, 1, 2, 2, 3, 3]
     ])
-    values_2d = xr.DataArray(values_val_2d)
+    agg = create_test_raster(data, backend)
+    return agg
 
-    values_val_3d = np.ones(4*3*8).reshape(3, 8, 4)
-    values_3d = xr.DataArray(
-        values_val_3d,
-        dims=['lat', 'lon', 'race']
-    )
-    values_3d['race'] = ['cat1', 'cat2', 'cat3', 'cat4']
 
+@pytest.fixture
+def data_values_3d(backend):
+    data = np.ones(4*3*8).reshape(3, 8, 4)
     if 'dask' in backend:
-        zones.data = da.from_array(zones.data, chunks=(3, 3))
-        values_2d.data = da.from_array(values_2d.data, chunks=(3, 3))
-        values_3d.data = da.from_array(values_3d.data, chunks=(3, 3, 1))
-
-    return zones, values_2d, values_3d
-
+        data = da.from_array(data, chunks=(3, 4, 2))
 
-def check_results(df_np, df_da, expected_results_dict):
-    # numpy case
-    assert isinstance(df_np, pd.DataFrame)
-    assert len(df_np.columns) == len(expected_results_dict)
+    agg = xr.DataArray(data, dims=['lat', 'lon', 'race'])
+    agg['race'] = ['cat1', 'cat2', 'cat3', 'cat4']
+    return agg
 
-    # zone column
-    assert (df_np['zone'] == expected_results_dict['zone']).all()
-
-    for col in df_np.columns[1:]:
-        assert np.isclose(
-            df_np[col], expected_results_dict[col], equal_nan=True
-        ).all()
-
-    if df_da is not None:
-        # dask case
-        assert isinstance(df_da, dd.DataFrame)
-        df_da = df_da.compute()
-        assert isinstance(df_da, pd.DataFrame)
-
-        # numpy results equal dask results, ignoring their indexes
-        assert np.array_equal(df_np.values, df_da.values, equal_nan=True)
 
-
-def test_stats():
-    # expected results
-    default_stats_results = {
+@pytest.fixture
+def result_default_stats():
+    expected_result = {
         'zone':  [0, 1, 2, 3],
         'mean':  [0, 1, 2, 2.4],
         'max':   [0, 1, 2, 3],
@@ -78,19 +60,13 @@ def test_stats():
         'var':   [0, 0, 0, 1.44],
         'count': [5, 6, 4, 5]
     }
+    return expected_result
 
-    # numpy case
-    zones_np, values_np, _ = create_zones_values(backend='numpy')
-    # default stats_funcs
-    df_np = stats(zones=zones_np, values=values_np)
-
-    # dask case
-    zones_da, values_da, _ = create_zones_values(backend='dask')
-    df_da = stats(zones=zones_da, values=values_da)
-    check_results(df_np, df_da, default_stats_results)
 
-    # expected results
-    stats_results_zone_0_3 = {
+@pytest.fixture
+def result_zone_ids_stats():
+    zone_ids = [0, 3]
+    expected_result = {
         'zone':  [0, 3],
         'mean':  [0, 2.4],
         'max':   [0, 3],
@@ -100,132 +76,166 @@ def test_stats():
         'var':   [0, 1.44],
         'count': [5, 5]
     }
+    return zone_ids, expected_result
 
-    # numpy case
-    df_np_zone_0_3 = stats(zones=zones_np, values=values_np, zone_ids=[0, 3])
 
-    # dask case
-    df_da_zone_0_3 = stats(zones=zones_da, values=values_da, zone_ids=[0, 3])
+def _double_sum(values):
+    return values.sum() * 2
 
-    check_results(df_np_zone_0_3, df_da_zone_0_3, stats_results_zone_0_3)
 
-    # ---- custom stats (NumPy only) ----
-    # expected results
-    custom_stats_results = {
+def _range(values):
+    return values.max() - values.min()
+
+
+@pytest.fixture
+def result_custom_stats():
+    zone_ids = [1, 2]
+    nodata_values = 0
+    expected_result = {
         'zone':       [1, 2],
         'double_sum': [12, 16],
         'range':      [0,   0],
     }
+    return nodata_values, zone_ids, expected_result
 
-    def _double_sum(values):
-        return values.sum() * 2
 
-    def _range(values):
-        return values.max() - values.min()
-
-    custom_stats = {
-        'double_sum': _double_sum,
-        'range': _range,
-    }
-
-    # numpy case
-    df_np = stats(
-        zones=zones_np, values=values_np, stats_funcs=custom_stats,
-        zone_ids=[1, 2], nodata_values=0
-    )
-    # dask case
-    df_da = None
-    check_results(df_np, df_da, custom_stats_results)
-
-
-def test_crosstab_2d():
-    # count agg, expected results
-    crosstab_2d_results = {
+@pytest.fixture
+def result_count_crosstab_2d():
+    zone_ids = [1, 2, 3]
+    cat_ids = [0, 1, 2]
+    expected_result = {
         'zone': [1, 2, 3],
         0:      [0, 0, 1],
         1:      [6, 0, 0],
         2:      [0, 4, 0],
     }
+    return zone_ids, cat_ids, expected_result
 
-    # numpy case
-    zones_np, values_np, _ = create_zones_values(backend='numpy')
-    df_np = crosstab(
-        zones=zones_np, values=values_np,
-        zone_ids=[1, 2, 3], cat_ids=[0, 1, 2],
-    )
-    # dask case
-    zones_da, values_da, _ = create_zones_values(backend='dask')
-    df_da = crosstab(
-        zones=zones_da, values=values_da, zone_ids=[1, 2, 3], nodata_values=3
-    )
-    check_results(df_np, df_da, crosstab_2d_results)
-
-    # percentage agg, expected results
 
-    crosstab_2d_percentage_results = {
+@pytest.fixture
+def result_percentage_crosstab_2d():
+    zone_ids = [1, 2]
+    cat_ids = [1, 2]
+    nodata_values = 3
+    expected_result = {
         'zone': [1,   2],
         1:      [100, 0],
         2:      [0,   100],
     }
-
-    # numpy case
-    df_np = crosstab(
-        zones=zones_np, values=values_np, zone_ids=[1, 2], cat_ids=[1, 2],
-        nodata_values=3, agg='percentage'
-    )
-    # dask case
-    df_da = crosstab(
-        zones=zones_da, values=values_da, zone_ids=[1, 2], cat_ids=[1, 2],
-        nodata_values=3, agg='percentage'
-    )
-    check_results(df_np, df_da, crosstab_2d_percentage_results)
+    return nodata_values, zone_ids, cat_ids, expected_result
 
 
-def test_crosstab_3d():
-    # expected results
-    crosstab_3d_results = {
+@pytest.fixture
+def result_crosstab_3d():
+    zone_ids = [1, 2, 3]
+    layer = -1
+    expected_result = {
         'zone': [1, 2, 3],
         'cat1': [6, 5, 6],
         'cat2': [6, 5, 6],
         'cat3': [6, 5, 6],
         'cat4': [6, 5, 6],
     }
+    return layer, zone_ids, expected_result
 
-    # numpy case
-    zones_np, _, values_np = create_zones_values(backend='numpy')
-    df_np = crosstab(
-        zones=zones_np, values=values_np, zone_ids=[1, 2, 3], layer=-1
-    )
-    # dask case
-    zones_da, _, values_da = create_zones_values(backend='dask')
-    df_da = crosstab(
-        zones=zones_da, values=values_da, zone_ids=[1, 2, 3],
-        cat_ids=['cat1', 'cat2', 'cat3', 'cat4'], layer=-1
-    )
-    check_results(df_np, df_da, crosstab_3d_results)
 
-    # ----- no values case ------
-    crosstab_3d_novalues_results = {
+@pytest.fixture
+def result_nodata_values_crosstab_3d():
+    zone_ids = [1, 2, 3]
+    layer = -1
+    nodata_values = 1
+    expected_result = {
         'zone': [1, 2, 3],
         'cat1': [0, 0, 0],
         'cat2': [0, 0, 0],
         'cat3': [0, 0, 0],
         'cat4': [0, 0, 0],
     }
+    return nodata_values, layer, zone_ids, expected_result
+
+
+def check_results(backend, df_result, expected_results_dict):
+    if 'dask' in backend:
+        # dask case, compute result
+        assert isinstance(df_result, dd.DataFrame)
+        df_result = df_result.compute()
+        assert isinstance(df_result, pd.DataFrame)
+
+    assert len(df_result.columns) == len(expected_results_dict)
+    # zone column
+    assert (df_result['zone'] == expected_results_dict['zone']).all()
+    # stats columns
+    for col in df_result.columns[1:]:
+        np.testing.assert_allclose(df_result[col], expected_results_dict[col])
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_default_stats(backend, data_zones, data_values_2d, result_default_stats):
+    df_result = stats(zones=data_zones, values=data_values_2d)
+    check_results(backend, df_result, result_default_stats)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_zone_ids_stats(backend, data_zones, data_values_2d, result_zone_ids_stats):
+    zone_ids, expected_result = result_zone_ids_stats
+    df_result = stats(zones=data_zones, values=data_values_2d, zone_ids=zone_ids)
+    check_results(backend, df_result, expected_result)
+
+
+@pytest.mark.parametrize("backend", ['numpy'])
+def test_custom_stats(backend, data_zones, data_values_2d, result_custom_stats):
+    # ---- custom stats (NumPy only) ----
+    custom_stats = {
+        'double_sum': _double_sum,
+        'range': _range,
+    }
+    nodata_values, zone_ids, expected_result = result_custom_stats
+    df_result = stats(
+        zones=data_zones, values=data_values_2d, stats_funcs=custom_stats,
+        zone_ids=zone_ids, nodata_values=nodata_values
+    )
+    check_results(backend, df_result, expected_result)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_count_crosstab_2d(backend, data_zones, data_values_2d, result_count_crosstab_2d):
+    zone_ids, cat_ids, expected_result = result_count_crosstab_2d
+    df_result = crosstab(
+        zones=data_zones, values=data_values_2d, zone_ids=zone_ids, cat_ids=cat_ids,
+    )
+    check_results(backend, df_result, expected_result)
+
 
-    # numpy case
-    zones_np, _, values_np = create_zones_values(backend='numpy')
-    df_np = crosstab(
-        zones=zones_np, values=values_np, layer=-1,
-        zone_ids=[1, 2, 3], nodata_values=1
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_percentage_crosstab_2d(backend, data_zones, data_values_2d, result_percentage_crosstab_2d):
+    nodata_values, zone_ids, cat_ids, expected_result = result_percentage_crosstab_2d
+    df_result = crosstab(
+        zones=data_zones, values=data_values_2d, zone_ids=zone_ids, cat_ids=cat_ids,
+        nodata_values=nodata_values, agg='percentage'
     )
-    # dask case
-    zones_da, _, values_da = create_zones_values(backend='dask')
-    df_da = crosstab(
-        zones=zones_da, values=values_da, layer=-1,
-        zone_ids=[1, 2, 3], nodata_values=1
+    check_results(backend, df_result, expected_result)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_crosstab_3d(backend, data_zones, data_values_3d, result_crosstab_3d):
+    layer, zone_ids, expected_result = result_crosstab_3d
+    df_result = crosstab(zones=data_zones, values=data_values_3d, zone_ids=zone_ids, layer=layer)
+    check_results(backend, df_result, expected_result)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_nodata_values_crosstab_3d(
+    backend,
+    data_zones,
+    data_values_3d,
+    result_nodata_values_crosstab_3d
+):
+    nodata_values, layer, zone_ids, expected_result = result_nodata_values_crosstab_3d
+    df_result = crosstab(
+        zones=data_zones, values=data_values_3d, zone_ids=zone_ids,
+        layer=layer, nodata_values=nodata_values
     )
-    check_results(df_np, df_da, crosstab_3d_novalues_results)
+    check_results(backend, df_result, expected_result)
 
 
 def test_apply():