diff --git a/xrspatial/aspect.py b/xrspatial/aspect.py
index abfe6b89..ff55e0f4 100644
--- a/xrspatial/aspect.py
+++ b/xrspatial/aspect.py
@@ -3,6 +3,8 @@
 import numpy as np
 import numba as nb
 
+import dask.array as da
+
 from numba import cuda
 
 from xarray import DataArray
@@ -147,7 +149,12 @@ def aspect(agg, name='aspect', use_cuda=True, pad=True, use_cupy=True):
         _horn_aspect_cuda[griddim, blockdim](data, out)
         if pad:
             out = out[pad_rows:-pad_rows, pad_cols:-pad_cols]
-        
+
+    elif isinstance(agg.data, da.Array):
+        out = agg.data.map_overlap(_horn_aspect,
+                                   depth=(1, 1),
+                                   boundary=np.nan,
+                                   meta=np.array(()))
     else:
         out = _horn_aspect(agg.data)
 
diff --git a/xrspatial/hillshade.py b/xrspatial/hillshade.py
index e9f047fd..2ffb45a3 100644
--- a/xrspatial/hillshade.py
+++ b/xrspatial/hillshade.py
@@ -1,9 +1,25 @@
-from __future__ import division, absolute_import
+from functools import partial
 
 import numpy as np
 
 from xarray import DataArray
 
+import dask.array as da
+
+
+def _hillshade(data, azimuth=225, angle_altitude=25):
+    azimuth = 360.0 - azimuth
+    x, y = np.gradient(data)
+    slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
+    aspect = np.arctan2(-x, y)
+    azimuthrad = azimuth*np.pi/180.
+    altituderad = angle_altitude*np.pi/180.
+    shaded = np.sin(altituderad) * np.sin(slope) + np.cos(altituderad) * np.cos(slope)*np.cos((azimuthrad - np.pi/2.) - aspect)
+    result = (shaded + 1) / 2
+    result[(0, -1), :] = np.nan
+    result[:, (0, -1)] = np.nan
+    return data
+
 
 def hillshade(agg, azimuth=225, angle_altitude=25, name='hillshade'):
     """Illuminates 2D DataArray from specific azimuth and altitude.
@@ -34,13 +50,15 @@ def hillshade(agg, azimuth=225, angle_altitude=25, name='hillshade'):
     Algorithm References:
      - http://geoexamples.blogspot.com/2014/03/shaded-relief-images-using-gdal-python.html
     """
-    azimuth = 360.0 - azimuth
-    x, y = np.gradient(agg.data)
-    slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
-    aspect = np.arctan2(-x, y)
-    azimuthrad = azimuth*np.pi/180.
-    altituderad = angle_altitude*np.pi/180.
-    shaded = np.sin(altituderad) * np.sin(slope) + np.cos(altituderad) * np.cos(slope)*np.cos((azimuthrad - np.pi/2.) - aspect)
-    data = (shaded + 1) / 2
-    return DataArray(data, name=name, dims=agg.dims,
+
+    if isinstance(agg.data, da.Array):
+        _func = partial(_hillshade, azimuth=azimuth, angle_altitude=angle_altitude)
+        out = agg.data.map_overlap(_func,
+                                   depth=(1, 1),
+                                   boundary=np.nan,
+                                   meta=np.array(()))
+    else:
+        out = _hillshade(agg.data, azimuth, angle_altitude)
+
+    return DataArray(out, name=name, dims=agg.dims,
                      coords=agg.coords, attrs=agg.attrs)
diff --git a/xrspatial/slope.py b/xrspatial/slope.py
index e8ce997f..342944f6 100644
--- a/xrspatial/slope.py
+++ b/xrspatial/slope.py
@@ -1,18 +1,32 @@
+# std lib
+from functools import partial
 from math import atan
-import numpy as np
-import numba as nb
+from typing import Union
+
+# 3rd-party
+try:
+    import cupy
+except ImportError:
+    class cupy(object):
+        ndarray = False
 
+import dask.array as da
+
+import numba as nb
 from numba import cuda
 
-from xarray import DataArray
+import numpy as np
+import xarray as xr
 
-from xrspatial.utils import ngjit
-from xrspatial.utils import has_cuda
+# local modules
 from xrspatial.utils import cuda_args
-
+from xrspatial.utils import get_dataarray_resolution
+from xrspatial.utils import has_cuda
+from xrspatial.utils import ngjit
+from xrspatial.utils import is_cupy_backed
 
 @ngjit
-def _horn_slope(data, cellsize_x, cellsize_y):
+def _cpu(data, cellsize_x, cellsize_y):
     out = np.zeros_like(data)
     out[:] = np.nan
     rows, cols = data.shape
@@ -33,8 +47,29 @@ def _horn_slope(data, cellsize_x, cellsize_y):
     return out
 
 
+def _run_numpy(data:np.ndarray,
+               cellsize_x:Union[int, float],
+               cellsize_y:Union[int, float]) -> np.ndarray:
+    out = _cpu(data, cellsize_x, cellsize_y)
+    return out
+
+
+def _run_dask_numpy(data:da.Array,
+                   cellsize_x:Union[int, float],
+                   cellsize_y:Union[int, float]) -> da.Array:
+    _func = partial(_cpu,
+                    cellsize_x=cellsize_x,
+                    cellsize_y=cellsize_y)
+
+    out = data.map_overlap(_func,
+                           depth=(1, 1),
+                           boundary=np.nan,
+                           meta=np.array(()))
+    return out
+
+
 @cuda.jit(device=True)
-def _gpu_slope(arr, cellsize_x, cellsize_y):
+def _gpu(arr, cellsize_x, cellsize_y):
     a = arr[2, 0]
     b = arr[2, 1]
     c = arr[2, 2]
@@ -53,28 +88,73 @@ def _gpu_slope(arr, cellsize_x, cellsize_y):
 
 
 @cuda.jit
-def _horn_slope_cuda(arr, cellsize_x_arr, cellsize_y_arr, out):
+def _run_gpu(arr, cellsize_x_arr, cellsize_y_arr, out):
     i, j = cuda.grid(2)
     di = 1
     dj = 1
     if (i-di >= 1 and i+di < out.shape[0] - 1 and 
         j-dj >= 1 and j+dj < out.shape[1] - 1):
-        out[i, j] = _gpu_slope(arr[i-di:i+di+1, j-dj:j+dj+1],
-                               cellsize_x_arr,
-                               cellsize_y_arr)
+        out[i, j] = _gpu(arr[i-di:i+di+1, j-dj:j+dj+1],
+                         cellsize_x_arr,
+                         cellsize_y_arr)
+
 
+def _run_cupy(data: cupy.ndarray,
+              cellsize_x: Union[int, float],
+              cellsize_y: Union[int, float]) -> cupy.ndarray:
 
-def slope(agg, name='slope', use_cuda=True, pad=True, use_cupy=True):
+    cellsize_x_arr = cupy.array([float(cellsize_x)], dtype='f4')
+    cellsize_y_arr = cupy.array([float(cellsize_y)], dtype='f4')
+
+    pad_rows = 3 // 2
+    pad_cols = 3 // 2
+    pad_width = ((pad_rows, pad_rows),
+                (pad_cols, pad_cols))
+
+    slope_data = np.pad(data, pad_width=pad_width, mode="reflect")
+
+    griddim, blockdim = cuda_args(slope_data.shape)
+    slope_agg = cupy.empty(slope_data.shape, dtype='f4')
+    slope_agg[:] = cupy.nan
+
+    _run_gpu[griddim, blockdim](slope_data,
+                                cellsize_x_arr,
+                                cellsize_y_arr,
+                                slope_agg)
+    out = slope_agg[pad_rows:-pad_rows, pad_cols:-pad_cols]
+    return out
+
+
+def _run_dask_cupy(data:da.Array,
+                   cellsize_x:Union[int, float],
+                   cellsize_y:Union[int, float]) -> da.Array:
+
+    msg = 'Upstream bug in dask prevents cupy backed arrays'
+    raise NotImplementedError(msg)
+
+    _func = partial(_run_cupy,
+                    cellsize_x=cellsize_x,
+                    cellsize_y=cellsize_y)
+
+    out = data.map_overlap(_func,
+                           depth=(1, 1),
+                           boundary=cupy.nan,
+                           dtype=cupy.float32,
+                           meta=cupy.array(()))
+    return out
+
+
+def slope(agg:xr.DataArray, name: str='slope') -> xr.DataArray:
     """Returns slope of input aggregate in degrees.
+
     Parameters
     ----------
-    agg : DataArray
+    agg : xr.DataArray
     name : str - name property of output xr.DataArray
-    use_cuda : bool - use CUDA device if available
 
     Returns
     -------
-    data: DataArray
+    data: xr.DataArray
 
     Notes:
     ------
@@ -85,61 +165,29 @@ def slope(agg, name='slope', use_cuda=True, pad=True, use_cupy=True):
       (Oxford University Press, New York), pp 406
     """
 
-    if not isinstance(agg, DataArray):
-        raise TypeError("agg must be instance of DataArray")
-
-    if not agg.attrs.get('res'):
-        raise ValueError('input xarray must have `res` attr.')
-
-    # get cellsize out from 'res' attribute
-    cellsize = agg.attrs.get('res')
-    if isinstance(cellsize, tuple) and len(cellsize) == 2 \
-            and isinstance(cellsize[0], (int, float)) \
-            and isinstance(cellsize[1], (int, float)):
-        cellsize_x, cellsize_y = cellsize
-    elif isinstance(cellsize, (int, float)):
-        cellsize_x = cellsize
-        cellsize_y = cellsize
-    else:
-        raise ValueError('`res` attr of input xarray must be a numeric'
-                         ' or a tuple of numeric values.')
+    cellsize_x, cellsize_y = get_dataarray_resolution(agg)
+
+    # numpy case
+    if isinstance(agg.data, np.ndarray):
+        out = _run_numpy(agg.data, cellsize_x, cellsize_y)
+
+    # cupy case
+    elif has_cuda() and isinstance(agg.data, cupy.ndarray):
+        out = _run_cupy(agg.data, cellsize_x, cellsize_y)
+
+    # dask + cupy case
+    elif has_cuda() and isinstance(agg.data, da.Array) and is_cupy_backed(agg):
+        out = _run_dask_cupy(agg.data, cellsize_x, cellsize_y)
     
-    if has_cuda() and use_cuda:
-        cellsize_x_arr = np.array([float(cellsize_x)], dtype='f4')
-        cellsize_y_arr = np.array([float(cellsize_y)], dtype='f4')
-
-        if pad:
-            pad_rows = 3 // 2
-            pad_cols = 3 // 2
-            pad_width = ((pad_rows, pad_rows),
-                        (pad_cols, pad_cols))
-        else:
-            # If padding is not desired, set pads to 0
-            pad_rows = 0
-            pad_cols = 0
-            pad_width = 0
-
-        slope_data = np.pad(agg.data, pad_width=pad_width, mode="reflect")
-
-        griddim, blockdim = cuda_args(slope_data.shape)
-        slope_agg = np.empty(slope_data.shape, dtype='f4')
-        slope_agg[:] = np.nan
-
-        if use_cupy:
-            import cupy
-            slope_agg = cupy.asarray(slope_agg)
-
-        _horn_slope_cuda[griddim, blockdim](slope_data,
-                                            cellsize_x_arr,
-                                            cellsize_y_arr,
-                                            slope_agg)
-        if pad:
-            slope_agg = slope_agg[pad_rows:-pad_rows, pad_cols:-pad_cols]
+    # dask + numpy case
+    elif isinstance(agg.data, da.Array):
+        out = _run_dask_numpy(agg.data, cellsize_x, cellsize_y)
+
     else:
-        slope_agg = _horn_slope(agg.data, cellsize_x, cellsize_y)
+        raise TypeError('Unsupported Array Type: {}'.format(type(agg.data)))
 
-    return DataArray(slope_agg,
-                     name=name,
-                     coords=agg.coords,
-                     dims=agg.dims,
-                     attrs=agg.attrs)
\ No newline at end of file
+    return xr.DataArray(out,
+                        name=name,
+                        coords=agg.coords,
+                        dims=agg.dims,
+                        attrs=agg.attrs)
\ No newline at end of file
diff --git a/xrspatial/tests/test_aspect.py b/xrspatial/tests/test_aspect.py
index de73ae98..b10daa88 100644
--- a/xrspatial/tests/test_aspect.py
+++ b/xrspatial/tests/test_aspect.py
@@ -1,6 +1,7 @@
-import xarray as xr
+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
@@ -325,4 +326,27 @@ def test_aspect_against_qgis_gpu():
     assert np.isclose(xrspatial_vals, qgis_vals, equal_nan=True).all()
 
     #assert (np.isnan(xrspatial_vals) | (
-    #        (0. <= xrspatial_vals) & (xrspatial_vals <= 360.))).all()
\ No newline at end of file
+    #        (0. <= xrspatial_vals) & (xrspatial_vals <= 360.))).all()
+
+
+def test_numpy_equals_dask():
+
+    # input data
+    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)
+
+    small_numpy_based_data_array = xr.DataArray(data, attrs={'res': (10.0, 10.0)})
+    small_das_based_data_array = xr.DataArray(da.from_array(data, chunks=(2, 2)),
+                                              attrs={'res': (10.0, 10.0)})
+
+    numpy_result = aspect(small_numpy_based_data_array, name='numpy_slope', use_cuda=False)
+    dask_result = aspect(small_das_based_data_array, name='dask_slope', use_cuda=False)
+    dask_result.data = dask_result.data.compute()
+
+    assert np.isclose(numpy_result, dask_result, equal_nan=True).all()
diff --git a/xrspatial/tests/test_hillshade.py b/xrspatial/tests/test_hillshade.py
index 1aba21bf..6f7360e4 100644
--- a/xrspatial/tests/test_hillshade.py
+++ b/xrspatial/tests/test_hillshade.py
@@ -3,6 +3,8 @@
 
 from xrspatial import hillshade
 
+import dask.array as da
+
 
 def _do_sparse_array(data_array):
     import random
@@ -29,6 +31,7 @@ def _do_gaussian_array():
 #
 # -----
 
+
 data_random = np.random.random_sample((100, 100))
 data_random_sparse = _do_sparse_array(data_random)
 data_gaussian = _do_gaussian_array()
@@ -47,3 +50,29 @@ def test_hillshade():
         assert np.all(da_gaussian_shade[coord] == da_gaussian[coord])
     assert da_gaussian_shade.mean() > 0
     assert da_gaussian_shade[60, 60] > 0
+
+
+def test_numpy_equals_dask():
+
+    # input data
+    data = 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)
+
+    attrs = {'res': (10.0, 10.0)}
+
+    small_numpy_based_data_array = xr.DataArray(data, attrs=attrs)
+    dask_data = da.from_array(data, chunks=(3, 3))
+    small_das_based_data_array = xr.DataArray(dask_data, attrs=attrs)
+
+    numpy_result = hillshade(small_numpy_based_data_array, name='numpy')
+    dask_result = hillshade(small_das_based_data_array, name='dask')
+    dask_result.data = dask_result.data.compute()
+
+    assert np.isclose(numpy_result, dask_result, equal_nan=True).all()
diff --git a/xrspatial/tests/test_slope.py b/xrspatial/tests/test_slope.py
index 380a9d31..901ec8da 100644
--- a/xrspatial/tests/test_slope.py
+++ b/xrspatial/tests/test_slope.py
@@ -1,130 +1,58 @@
-
 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
 
 
-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():
-    _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_random = np.random.random_sample((100, 100))
-data_random_sparse = _do_sparse_array(data_random)
-data_gaussian = _do_gaussian_array()
-
-
-def test_invalid_res_attr():
-    """
-    Assert 'res' attribute of input xarray
-    """
-    da = xr.DataArray(data_gaussian, attrs={})
-    with pytest.raises(ValueError) as e_info:
-        da_slope = slope(da)
-        assert e_info
-
-    da = xr.DataArray(data_gaussian, attrs={'res': 'any_string'})
-    with pytest.raises(ValueError) as e_info:
-        da_slope = slope(da)
-        assert e_info
-
-    da = xr.DataArray(data_gaussian, attrs={'res': ('str_tuple', 'str_tuple')})
-    with pytest.raises(ValueError) as e_info:
-        da_slope = slope(da)
-        assert e_info
-
-    da = xr.DataArray(data_gaussian, attrs={'res': (1, 2, 3)})
-    with pytest.raises(ValueError) as e_info:
-        da_slope = slope(da)
-        assert e_info
-
-
-def test_slope_transfer_function():
-    """
-    Assert slope transfer function
-    """
-    da = xr.DataArray(data_gaussian, attrs={'res': 1})
-    da_slope = slope(da)
-    # default name
-    assert da_slope.name == 'slope'
-    assert da_slope.dims == da.dims
-    assert da_slope.attrs == da.attrs
-    assert da.shape == da_slope.shape
-    for coord in da.coords:
-        assert np.all(da_slope[coord] == da[coord])
-
-    assert da_slope.sum() > 0
-
-    # In the middle of the array, there is the maximum of the gaussian;
-    # And there the slope must be zero.
-    _imax = np.where(da == da.max())
-    assert da_slope[_imax] == 0
-
-    # same result when cellsize_x = cellsize_y = 1
-    da = xr.DataArray(data_gaussian, attrs={'res': (1.0, 1.0)})
-    da_slope = slope(da)
-    assert da_slope.dims == da.dims
-    assert da_slope.attrs == da.attrs
-    assert da.shape == da_slope.shape
-    for coord in da.coords:
-        assert np.all(da_slope[coord] == da[coord])
-
-    assert da_slope.sum() > 0
-
-    # In the middle of the array, there is the maximum of the gaussian;
-    # And there the slope must be zero.
-    _imax = np.where(da == da.max())
-    assert da_slope[_imax] == 0
-
+# 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)
+
+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)
+
+elevation2 = 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_slope_against_qgis():
-    # input data
-    data = 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)
-    small_da = xr.DataArray(data, attrs={'res': (10.0, 10.0)})
 
-    # slope by QGIS
-    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)
+    small_da = xr.DataArray(elevation, attrs={'res': (10.0, 10.0)})
 
     # slope by xrspatial
-    xrspatial_slope = slope(small_da, name='slope_agg', use_cuda=False)
+    xrspatial_slope = slope(small_da, name='slope_agg')
 
     # validate output attributes
     assert xrspatial_slope.dims == small_da.dims
@@ -144,37 +72,17 @@ def test_slope_against_qgis():
 
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
 def test_slope_against_qgis_gpu():
-    # input data
-    data = 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)
-    small_da = xr.DataArray(data, attrs={'res': (10.0, 10.0)})
 
-    # slope by QGIS
-    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)
+    import cupy
 
-    # slope by xrspatial
-    xrspatial_slope = slope(small_da, name='slope_agg', use_cuda=True)
+    small_da = xr.DataArray(elevation, attrs={'res': (10.0, 10.0)})
+    small_da_cupy = xr.DataArray(cupy.asarray(elevation), attrs={'res': (10.0, 10.0)})
+    xrspatial_slope = slope(small_da_cupy, name='slope_cupy')
 
     # validate output attributes
     assert xrspatial_slope.dims == small_da.dims
     assert xrspatial_slope.attrs == small_da.attrs
     assert xrspatial_slope.shape == small_da.shape
-    assert xrspatial_slope.name == 'slope_agg'
     for coord in small_da.coords:
         assert np.all(xrspatial_slope[coord] == small_da[coord])
 
@@ -189,20 +97,76 @@ def test_slope_against_qgis_gpu():
 @pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
 def test_slope_gpu_equals_cpu():
 
-    # input data
-    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)
+    import cupy
+
+    small_da = xr.DataArray(elevation2, attrs={'res': (10.0, 10.0)})
+    cpu = slope(small_da, name='numpy_result')
+
+    small_da_cupy = xr.DataArray(cupy.asarray(elevation2), attrs={'res': (10.0, 10.0)})
+    gpu = slope(small_da_cupy, name='cupy_result')
+    assert isinstance(gpu.data, cupy.ndarray)
+
+    assert np.isclose(cpu, gpu, equal_nan=True).all()
+
+
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+def _dask_cupy_equals_numpy_cpu():
+
+    # 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
+
+    import cupy
+
+    cupy_data = cupy.asarray(elevation2)
+    dask_cupy_data = da.from_array(cupy_data, chunks=(3, 3))
+
+    small_da = xr.DataArray(elevation2, attrs={'res': (10.0, 10.0)})
+    cpu = slope(small_da, name='numpy_result')
+
+    small_dask_cupy = xr.DataArray(dask_cupy_data, attrs={'res': (10.0, 10.0)})
+    gpu = slope(small_dask_cupy, name='cupy_result')
+
+    assert np.isclose(cpu, gpu, equal_nan=True).all()
+
+
+def test_slope_numpy_equals_dask():
+
+    small_numpy_based_data_array = xr.DataArray(elevation2, attrs={'res': (10.0, 10.0)})
+    small_das_based_data_array = xr.DataArray(da.from_array(elevation2, chunks=(3, 3)),
+                                              attrs={'res': (10.0, 10.0)})
 
+    numpy_slope = slope(small_numpy_based_data_array, name='numpy_slope')
+    dask_slope = slope(small_das_based_data_array, name='dask_slope')
+    assert isinstance(dask_slope.data, da.Array)
+
+    dask_slope.data = dask_slope.data.compute()
+
+    assert np.isclose(numpy_slope, dask_slope, equal_nan=True).all()
+
+
+def test_slope_with_dask_array():
+
+    import dask.array as da
+
+    data = da.from_array(elevation, chunks=(3, 3))
     small_da = xr.DataArray(data, attrs={'res': (10.0, 10.0)})
 
-    # aspect by xrspatial
-    cpu = slope(small_da, name='aspect_agg', use_cuda=False)
-    gpu = slope(small_da, name='aspect_agg', use_cuda=True)
+    # slope by xrspatial
+    xrspatial_slope = slope(small_da, name='slope_agg')
+    xrspatial_slope.data = xrspatial_slope.data.compute()
 
-    assert np.isclose(cpu, gpu, equal_nan=True).all()
\ No newline at end of file
+    # validate output attributes
+    assert xrspatial_slope.dims == small_da.dims
+    assert xrspatial_slope.attrs == small_da.attrs
+    assert xrspatial_slope.shape == small_da.shape
+    assert xrspatial_slope.name == 'slope_agg'
+    for coord in small_da.coords:
+        assert np.all(xrspatial_slope[coord] == small_da[coord])
+
+    # validate output values
+    # ignore border edges
+    xrspatial_vals = xrspatial_slope.values[1:-1, 1:-1]
+    qgis_vals = qgis_slope[1:-1, 1:-1]
+    assert (np.isclose(xrspatial_vals, qgis_vals, equal_nan=True).all() | (
+                np.isnan(xrspatial_vals) & np.isnan(qgis_vals))).all()
diff --git a/xrspatial/utils.py b/xrspatial/utils.py
index ebe39ee4..ce69be63 100644
--- a/xrspatial/utils.py
+++ b/xrspatial/utils.py
@@ -1,6 +1,7 @@
 from math import ceil
 import numba as nb
 import numpy as np
+import xarray as xr
 
 from numba import cuda
 
@@ -61,6 +62,48 @@ def cuda_args(shape):
     return bpg, tpb
 
 
+def is_cupy_backed(agg: xr.DataArray):
+    return type(agg.data._meta).__module__.split('.')[0] == 'cupy'
+
+
+def calc_res(raster):
+    """Calculate the resolution of xarray.DataArray raster and return it as the
+    two-tuple (xres, yres).
+
+    Notes
+    -----
+
+    Sourced from datashader.utils
+    """
+    h, w = raster.shape[-2:]
+    ydim, xdim = raster.dims[-2:]
+    xcoords = raster[xdim].values
+    ycoords = raster[ydim].values
+    xres = (xcoords[-1] - xcoords[0]) / (w - 1)
+    yres = (ycoords[0] - ycoords[-1]) / (h - 1)
+    return xres, yres
+
+
+def get_dataarray_resolution(agg: xr.DataArray):
+
+    if not agg.attrs.get('res'):
+        raise ValueError('input xarray must have `res` attr.')
+
+    # get cellsize out from 'res' attribute
+    cellsize = agg.attrs.get('res')
+    if isinstance(cellsize, tuple) and len(cellsize) == 2 \
+            and isinstance(cellsize[0], (int, float)) \
+            and isinstance(cellsize[1], (int, float)):
+        cellsize_x, cellsize_y = cellsize
+    elif isinstance(cellsize, (int, float)):
+        cellsize_x = cellsize
+        cellsize_y = cellsize
+    else:
+        cellsize_x, cellsize_y = calc_res(agg)
+    
+    return cellsize_x, cellsize_y
+
+
 def lnglat_to_meters(longitude, latitude):
     """
     Projects the given (longitude, latitude) values into Web Mercator