Added focal_stats (#453)

* commented out warnings

* removed calc_mean and calc_sum from documentation

* added focal_stats

* added stats functions

* added tests

* corrected docstrings

* flake8 fixes

* convolution_2d

* flake8 fixes

* remove commented code

* fix focal docstrings

* add convolution_2d and focal_stats to docs

Co-authored-by: giancastro <giancastrok@gmail.com>

Author: thuydotm

docs/source/reference/focal.rst

@@ -31,9 +31,9 @@ Focal Statistics
 .. autosummary::
    :toctree: _autosummary
 
+   xrspatial.convolution.convolution_2d
    xrspatial.convolution.annulus_kernel
    xrspatial.convolution.calc_cellsize
-   xrspatial.focal.calc_mean
-   xrspatial.focal.calc_sum
    xrspatial.convolution.circle_kernel
    xrspatial.focal.custom_kernel
+   xrspatial.focal.focal_stats

xrspatial/convolution.py

@@ -1,10 +1,10 @@
 from functools import partial
 
 import re
-import warnings
 
 import numpy as np
 import dask.array as da
+from xarray import DataArray
 
 from numba import cuda, float32, prange, jit
 
@@ -52,13 +52,8 @@ def _get_distance(distance_str):
         raise ValueError("Invalid distance.")
 
     unit = DEFAULT_UNIT
-    if len(splits) == 1:
-        with warnings.catch_warnings():
-            warnings.simplefilter('default')
-            warnings.warn('Raster distance unit not provided. '
-                          'Use meter as default.', Warning)
 
-    elif len(splits) == 2:
+    if len(splits) == 2:
         unit = splits[1]
 
     number = splits[0]
@@ -163,10 +158,6 @@ def calc_cellsize(raster):
         unit = raster.attrs['unit']
     else:
         unit = DEFAULT_UNIT
-        with warnings.catch_warnings():
-            warnings.simplefilter('default')
-            warnings.warn('Raster distance unit not provided. '
-                          'Use meter as default.', Warning)
 
     cellsize_x, cellsize_y = get_dataarray_resolution(raster)
     cellsize_x = _to_meters(cellsize_x, unit)
@@ -286,25 +277,6 @@ def annulus_kernel(cellsize_x, cellsize_y, outer_radius, inner_radius):
          [0., 1., 1., 1., 1., 0., 1., 1., 1., 1., 0.],
          [0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.]])
     """
-    # validate radii, convert to meters
-    r2 = _get_distance(str(outer_radius))
-    r1 = _get_distance(str(inner_radius))
-
-    # Validate that outer radius is indeed outer radius
-    if r2 > r1:
-        r_outer = r2
-        r_inner = r1
-    else:
-        r_outer = r1
-        r_inner = r2
-
-    if r_outer - r_inner < np.sqrt((cellsize_x / 2)**2 + (cellsize_y / 2)**2):
-        with warnings.catch_warnings():
-            warnings.simplefilter('default')
-            warnings.warn(
-                'Annulus radii are closer than cellsize distance.', Warning
-            )
-
     # Get the two circular kernels for the annulus
     kernel_outer = circle_kernel(cellsize_x, cellsize_y, outer_radius)
     kernel_inner = circle_kernel(cellsize_x, cellsize_y, inner_radius)
@@ -534,3 +506,35 @@ def convolve_2d(data, kernel):
         raise TypeError('Unsupported Array Type: {}'.format(type(data)))
 
     return out
+
+
+def convolution_2d(agg, kernel):
+    """
+    Calculates, for all inner cells of an array, the 2D convolution of
+    each cell via Numba. To account for edge cells, a pad can be added
+    to the image array. Convolution is frequently used for image
+    processing, such as smoothing, sharpening, and edge detection of
+    images by eliminating spurious data or enhancing features in the
+    data.
+
+    Parameters
+    ----------
+    agg : xarray.DataArray
+        2D array of values to processed and padded.
+    kernel : array-like object
+        Impulse kernel, determines area to apply impulse function for
+        each cell.
+
+    Returns
+    -------
+    convolve_agg : xarray.DataArray
+        2D array representation of the impulse function.
+    """
+
+    # wrapper of convolve_2d
+    out = convolve_2d(agg.data, kernel)
+
+    return DataArray(out,
+                     coords=agg.coords,
+                     dims=agg.dims,
+                     attrs=agg.attrs)

xrspatial/focal.py

@@ -1,19 +1,19 @@
-from functools import partial
+import numpy as np
+import pandas as pd
+import xarray as xr
+import dask.array as da
 
+from functools import partial
 from math import isnan
-
-from numba import prange
-import numpy as np
+from numba import prange, cuda
 from xarray import DataArray
-import dask.array as da
 
 try:
     import cupy
 except ImportError:
     class cupy(object):
         ndarray = False
 
-from numba import cuda
 from xrspatial.utils import cuda_args
 from xrspatial.utils import has_cuda
 from xrspatial.utils import is_cupy_backed
@@ -208,119 +208,40 @@ def mean(agg, passes=1, excludes=[np.nan], name='mean'):
 
 
 @ngjit
-def calc_mean(array):
-    """
-    Calculates the mean of an array.
+def _calc_mean(array):
+    return np.nanmean(array)
 
-    Parameters
-    ----------
-    array : numpy.Array
-        Array of input values.
 
-    Returns
-    -------
-    array_sum : float
-        Mean of input data.
+@ngjit
+def _calc_sum(array):
+    return np.nansum(array)
 
-    Examples
-    --------
-    .. sourcecode:: python
 
-        >>> from xrspatial.focal import calc_mean
-        >>> import numpy as np
-
-        >>> # 1D Array of Integers
-        >>> array1 = np.array([1, 2, 3, 4, 5])
-        >>> print(array1)
-        [1 2 3 4 5]
-
-        >>> # Calculate Mean
-        >>> array_mean = calc_mean(array1)
-        >>> print(array_mean)
-        3.0
-
-        >>> # 2D Array of Floats
-        >>> array2 = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
-        >>> print(array2)
-        [[1. 2. 3.]
-         [4. 5. 6.]]
-
-        >>> # Calculate Mean
-        >>> array_mean = calc_mean(array2)
-        >>> print(array_mean)
-        3.5
-
-        >>> # 3D Array of Integers
-        >>> array3 = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
-        >>> print(array3)
-        [[1 2 3]
-         [4 5 6]
-         [7 8 9]]
-
-        >>> # Calculate Mean
-        >>> array_mean = calc_mean(array3)
-        >>> print(array_mean)
-        5.0
-    """
-    return np.nanmean(array)
+@ngjit
+def _calc_min(array):
+    return np.nanmin(array)
 
 
 @ngjit
-def calc_sum(array):
-    """
-    Calculates the sum of an array.
+def _calc_max(array):
+    return np.nanmax(array)
 
-    Parameters
-    ----------
-    array : numpy.Array
-        Array of input values.
 
-    Returns
-    -------
-    array_sum : float
-        Sum of input data.
+@ngjit
+def _calc_std(array):
+    return np.nanstd(array)
 
-    Examples
-    --------
-    .. sourcecode:: python
 
-        >>> from xrspatial.focal import calc_sum
-        >>> import numpy as np
-
-        >>> # 1D Array of Integers
-        >>> array1 = np.array([1, 2, 3, 4, 5])
-        >>> print(array1)
-        [1 2 3 4 5]
-
-        >>> # Calculate Sum
-        >>> array_sum = calc_sum(array1)
-        >>> print(array_sum)
-        15
-
-        >>> # 2D Array of Floats
-        >>> array2 = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
-        >>> print(array2)
-        [[1. 2. 3.]
-         [4. 5. 6.]]
-
-        >>> # Calculate Sum
-        >>> array_sum = calc_sum(array2)
-        >>> print(array_sum)
-        21.0
-
-        >>> # 3D Array of Integers
-        >>> array3 = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
-        >>> print(array3)
-        [[1 2 3]
-         [4 5 6]
-         [7 8 9]]
-
-        >>> # Calculate Sum
-        >>> agg_sum = calc_sum(array3)
-        >>> print(agg_sum)
-        45
-    """
-    return np.nansum(array)
+@ngjit
+def _calc_range(array):
+    value_min = _calc_min(array)
+    value_max = _calc_max(array)
+    return value_max - value_min
+
+
+@ngjit
+def _calc_var(array):
+    return np.nanvar(array)
 
 
 @ngjit
@@ -419,17 +340,17 @@ def _apply(data, kernel, func):
     return out
 
 
-def apply(raster, kernel, func=calc_mean):
+def apply(raster, kernel, func=_calc_mean):
     """
-    Returns Mean filtered array using a user-created window.
+    Returns custom function applied array using a user-created window.
 
     Parameters
     ----------
     raster : xarray.DataArray
         2D array of input values to be filtered.
-    kernel : Numpy Array
+    kernel : numpy.array
         2D array where values of 1 indicate the kernel.
-    func : xrspatial.focal.calc_mean
+    func : callable, default=xrspatial.focal._calc_mean
         Function which takes an input array and returns an array.
 
     Returns
@@ -547,6 +468,52 @@ def apply(raster, kernel, func=calc_mean):
     return result
 
 
+def focal_stats(agg,
+                kernel,
+                stats_funcs=[
+                    'mean', 'max', 'min', 'range', 'std', 'var', 'sum'
+                ]):
+    """
+    Calculates statistics of the values within a specified focal neighborhood
+    for each pixel in an input raster. The statistics types are Mean, Maximum,
+    Minimum, Range, Standard deviation, Variation and Sum.
+
+    Parameters
+    ----------
+    agg : xarray.DataArray
+        2D array of input values to be analysed.
+    kernel : numpy.array
+        2D array where values of 1 indicate the kernel.
+    stats_funcs: list of string
+        List of statistics types to be calculated.
+        Default set to ['mean', 'max', 'min', 'range', 'std', 'var', 'sum'].
+
+    Returns
+    -------
+    stats_agg : xarray.DataArray of same type as `agg`
+        3D array with dimensions of `(stat, y, x)` and with values
+        indicating the focal stats.
+    """
+
+    _function_mapping = {
+        'mean': _calc_mean,
+        'max': _calc_max,
+        'min': _calc_min,
+        'range': _calc_range,
+        'std': _calc_std,
+        'var': _calc_var,
+        'sum': _calc_sum
+    }
+
+    stats_aggs = []
+    for stats in stats_funcs:
+        stats_agg = apply(agg, kernel, func=_function_mapping[stats])
+        stats_aggs.append(stats_agg)
+
+    stats = xr.concat(stats_aggs, pd.Index(stats_funcs, name='stats'))
+    return stats
+
+
 @ngjit
 def _calc_hotspots_numpy(z_array):
     out = np.zeros_like(z_array, dtype=np.int8)

xrspatial/tests/test_focal.py

@@ -8,7 +8,7 @@
 from xrspatial.utils import ngjit
 
 from xrspatial import mean
-from xrspatial.focal import hotspots, apply
+from xrspatial.focal import hotspots, apply, focal_stats
 from xrspatial.convolution import (
     convolve_2d, calc_cellsize, circle_kernel, annulus_kernel
 )
@@ -287,6 +287,27 @@ def func_zero_cpu(x):
         assert e_info
 
 
+def test_focal_stats_cpu():
+    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)
+
+    numpy_focalstats = focal_stats(numpy_agg, kernel)
+    assert isinstance(numpy_focalstats.data, np.ndarray)
+    assert numpy_focalstats.ndim == 3
+    assert numpy_agg.shape == numpy_focalstats.shape[1:]
+
+    dask_numpy_focalstats = focal_stats(dask_numpy_agg, kernel)
+    assert isinstance(dask_numpy_focalstats.data, da.Array)
+
+    assert np.isclose(
+        numpy_focalstats, dask_numpy_focalstats.compute(), equal_nan=True
+    ).all()
+
+
 def test_hotspot():
     n, m = 10, 10
     data = np.zeros((n, m), dtype=float)