fixed UNKNOWN CUDA error on slope (#141)

* fixed UNKNOWN CUDA error on slope

* added cupy support

* aspect wip

* gpu aspect tests passing

* added gpu == cpu slope test

* added skipif CUDA Device not Available

* added GPU support for SAVI

* fixed aspect tests

Author: brendancol

xrspatial/aspect.py

@@ -1,15 +1,23 @@
+from math import atan2
+from math import pi
 import numpy as np
+import numba as nb
+
+from numba import cuda
+
 from xarray import DataArray
 
 from xrspatial.utils import ngjit
-
+from xrspatial.utils import has_cuda
+from xrspatial.utils import cuda_args
 
 RADIAN = 180 / np.pi
 
 
 @ngjit
 def _horn_aspect(data):
     out = np.zeros_like(data, dtype=np.float64)
+    out[:] = np.nan
     rows, cols = data.shape
     for y in range(1, rows-1):
         for x in range(1, cols-1):
@@ -28,7 +36,7 @@ def _horn_aspect(data):
 
             if dz_dx == 0 and dz_dy == 0:
                 # flat surface, slope = 0, thus invalid aspect
-                out[y, x] = np.nan
+                out[y, x] = -1.
             else:
                 aspect = np.arctan2(dz_dy, -dz_dx) * RADIAN
                 # convert to compass direction values (0-360 degrees)
@@ -42,7 +50,57 @@ def _horn_aspect(data):
     return out
 
 
-def aspect(agg, name='aspect'):
+@cuda.jit(device=True)
+def _gpu_aspect(arr):
+
+    a = arr[0, 0]
+    b = arr[0, 1]
+    c = arr[0, 2]
+    d = arr[1, 0]
+    f = arr[1, 2]
+    g = arr[2, 0]
+    h = arr[2, 1]
+    i = arr[2, 2]
+
+    two = nb.int32(2.)  # reducing size to int8 causes wrong results
+    eight = nb.int32(8.)  # reducing size to int8 causes wrong results
+    ninety = nb.float32(90.)
+
+    dz_dx = ((c + two * f + i) - (a + two * d + g)) / eight
+    dz_dy = ((g + two * h + i) - (a + two * b + c)) / eight
+
+    if dz_dx == 0 and dz_dy == 0:
+        # flat surface, slope = 0, thus invalid aspect
+        aspect = nb.float32(-1.)  # TODO: return null instead
+    else:
+        aspect = atan2(dz_dy, -dz_dx) * nb.float32(57.29578)
+        # convert to compass direction values (0-360 degrees)
+        if aspect < nb.float32(0.):
+            aspect = ninety - aspect
+        elif aspect > ninety:
+            aspect = nb.float32(360.0) - aspect + ninety
+        else:
+            aspect = ninety - aspect
+    
+    if aspect > nb.float32(359.999):  # lame float equality check...
+        return nb.float32(0.)
+    else:
+        return aspect
+
+
+@cuda.jit
+def _horn_aspect_cuda(arr, out):
+    minus_one = nb.float32(-1.)
+    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_aspect(arr[i-di:i+di+1, j-dj:j+dj+1])
+
+
+
+def aspect(agg, name='aspect', use_cuda=True, pad=True, use_cupy=True):
     """Returns downward slope direction in compass degrees (0 - 360) with 0 at 12 o'clock.
 
     Parameters
@@ -63,7 +121,37 @@ def aspect(agg, name='aspect'):
     if not isinstance(agg, DataArray):
         raise TypeError("agg must be instance of DataArray")
 
-    return DataArray(_horn_aspect(agg.data),
+    if has_cuda() and use_cuda:
+
+        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
+
+        data = np.pad(agg.data, pad_width=pad_width, mode="reflect")
+
+        griddim, blockdim = cuda_args(data.shape)
+        out = np.empty(data.shape, dtype='f4')
+        out[:] = np.nan
+
+        if use_cupy:
+            import cupy
+            out = cupy.asarray(out)
+
+        _horn_aspect_cuda[griddim, blockdim](data, out)
+        if pad:
+            out = out[pad_rows:-pad_rows, pad_cols:-pad_cols]
+        
+    else:
+        out = _horn_aspect(agg.data)
+
+    return DataArray(out,
                      name=name,
                      dims=agg.dims,
                      coords=agg.coords,

xrspatial/multispectral.py

@@ -1,11 +1,16 @@
 import numpy as np
+import numba as nb
+
+from numba import cuda
 
 import datashader as ds
 
 from PIL import Image
 
 from xarray import DataArray
 
+from xrspatial.utils import has_cuda
+from xrspatial.utils import cuda_args
 from xrspatial.utils import ngjit
 
 
@@ -408,7 +413,19 @@ def _savi(nir_data, red_data, soil_factor):
     return out
 
 
-def savi(nir_agg, red_agg, soil_factor=1.0, name='savi'):
+@cuda.jit
+def _savi_gpu(nir_data, red_data, soil_factor, out):
+    y, x = cuda.grid(2)
+    if y < out.shape[0] and x < out.shape[1]:
+        nir = nir_data[y, x]
+        red = red_data[y, x]
+        numerator = nir - red
+        soma = nir + red + soil_factor
+        denominator = soma * (1.0 + soil_factor)
+        out[y, x] = numerator / denominator
+
+
+def savi(nir_agg, red_agg, soil_factor=1.0, name='savi', use_cuda=True, use_cupy=True):
     """Returns Soil Adjusted Vegetation Index (SAVI).
 
     Parameters
@@ -441,7 +458,26 @@ def savi(nir_agg, red_agg, soil_factor=1.0, name='savi'):
     if soil_factor > 1.0 or soil_factor < -1.0:
         raise ValueError("soil factor must be between (-1.0, 1.0)")
 
-    return DataArray(_savi(nir_agg.data, red_agg.data, soil_factor),
+    nir_data = nir_agg.data
+    red_data = red_agg.data
+
+    if has_cuda() and use_cuda:
+        griddim, blockdim = cuda_args(nir_data.shape)
+        out = np.empty(nir_data.shape, dtype='f4')
+        out[:] = np.nan
+
+        if use_cupy:
+            import cupy
+            out = cupy.asarray(out)
+
+        _savi_gpu[griddim, blockdim](nir_data,
+                                     red_data,
+                                     soil_factor,
+                                     out)
+    else:
+        out = _savi(nir_agg.data, red_agg.data, soil_factor)
+
+    return DataArray(out,
                      name=name,
                      coords=nir_agg.coords,
                      dims=nir_agg.dims,

xrspatial/slope.py

@@ -1,5 +1,6 @@
 from math import atan
 import numpy as np
+import numba as nb
 
 from numba import cuda
 
@@ -13,6 +14,7 @@
 @ngjit
 def _horn_slope(data, cellsize_x, cellsize_y):
     out = np.zeros_like(data)
+    out[:] = np.nan
     rows, cols = data.shape
     for y in range(1, rows-1):
         for x in range(1, cols-1):
@@ -42,27 +44,27 @@ def _gpu_slope(arr, cellsize_x, cellsize_y):
     h = arr[0, 1]
     i = arr[0, 2]
 
-    dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cellsize_x[0])
-    dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cellsize_y[0])
-    p = (dz_dx * dz_dx + dz_dy * dz_dy) ** .5
-    return atan(p) * 57.29578
+    two = nb.int32(2.)  # reducing size to int8 causes wrong results
+
+    dz_dx = ((c + two * f + i) - (a + two * d + g)) / (nb.float32(8.) * cellsize_x[0])
+    dz_dy = ((g + two * h + i) - (a + two * b + c)) / (nb.float32(8.) * cellsize_y[0])
+    p = (dz_dx * dz_dx + dz_dy * dz_dy) ** nb.float32(.5)
+    return atan(p) * nb.float32(57.29578)
 
 
 @cuda.jit
 def _horn_slope_cuda(arr, cellsize_x_arr, cellsize_y_arr, out):
     i, j = cuda.grid(2)
     di = 1
     dj = 1
-    if (i-di >= 0 and i+di < out.shape[0] and 
-        j-dj >= 0 and j+dj < out.shape[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)
-    else:
-        out[i, j] = np.nan
 
 
-def slope(agg, name='slope', use_cuda=True):
+def slope(agg, name='slope', use_cuda=True, pad=True, use_cupy=True):
     """Returns slope of input aggregate in degrees.
     Parameters
     ----------
@@ -103,21 +105,41 @@ def slope(agg, name='slope', use_cuda=True):
                          ' or a tuple of numeric values.')
 
     if has_cuda() and use_cuda:
-        cellsize_x_arr =  np.array([float(cellsize_x)], dtype='f8')
-        cellsize_y_arr =  np.array([float(cellsize_y)], dtype='f8')
-
-        griddim, blockdim = cuda_args(agg.data.shape)
-        slope_agg = np.zeros(agg.data.shape, dtype='f8')
-        _horn_slope_cuda[griddim, blockdim](agg.data,
+        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)
-        pass
+        if pad:
+            slope_agg = slope_agg[pad_rows:-pad_rows, pad_cols:-pad_cols]
     else:
         slope_agg = _horn_slope(agg.data, cellsize_x, cellsize_y)
 
     return DataArray(slope_agg,
                      name=name,
                      coords=agg.coords,
                      dims=agg.dims,
-                     attrs=agg.attrs)
+                     attrs=agg.attrs)