import copy
from functools import partial
from math import isnan, sqrt
import dask.array as da
import numba as nb
import numpy as np
import pandas as pd
import xarray as xr
from numba import cuda, prange
from xarray import DataArray
try:
import cupy
except ImportError:
class cupy(object):
ndarray = False
from xrspatial.convolution import convolve_2d, custom_kernel
from xrspatial.utils import ArrayTypeFunctionMapping, cuda_args, ngjit, not_implemented_func
# TODO: Make convolution more generic with numba first-class functions.
@ngjit
def _equal_numpy(x, y):
if x == y or (np.isnan(x) and np.isnan(y)):
return True
return False
@ngjit
def _mean_numpy(data, excludes):
out = np.zeros_like(data)
rows, cols = data.shape
for y in range(rows):
for x in range(cols):
exclude = False
for ex in excludes:
if _equal_numpy(data[y, x], ex):
exclude = True
break
if not exclude:
left = max(x-1, 0)
right = min(x+2, cols)
bottom = max(y-1, 0)
top = min(y+2, rows)
kernel_data = data[bottom:top, left:right]
out[y, x] = np.nanmean(kernel_data)
else:
out[y, x] = data[y, x]
return out
def _mean_dask_numpy(data, excludes):
_func = partial(_mean_numpy, excludes=excludes)
out = data.map_overlap(_func,
depth=(1, 1),
boundary=np.nan,
meta=np.array(()))
return out
@cuda.jit
def _mean_gpu(data, excludes, out):
i, j = cuda.grid(2)
for ex in excludes:
if (data[i, j] == ex) or (isnan(data[i, j]) and isnan(ex)):
out[i, j] = data[i, j]
return
rows, cols = out.shape
if 0 <= i < rows and 0 <= j < cols:
left = max(j - 1, 0)
right = min(j + 2, cols)
bottom = max(i - 1, 0)
top = min(i + 2, rows)
sum = 0
num = 0
for y in range(bottom, top):
for x in range(left, right):
if not isnan(data[y, x]):
sum += data[y, x]
num += 1
if num > 0:
out[i, j] = sum / num
def _mean_cupy(data, excludes):
griddim, blockdim = cuda_args(data.shape)
out = cupy.empty(data.shape, dtype='f4')
out[:] = cupy.nan
_mean_gpu[griddim, blockdim](data, cupy.asarray(excludes), out)
return out
def _mean(data, excludes):
agg = xr.DataArray(data)
mapper = ArrayTypeFunctionMapping(
numpy_func=_mean_numpy,
cupy_func=_mean_cupy,
dask_func=_mean_dask_numpy,
dask_cupy_func=lambda *args: not_implemented_func(
*args, messages='mean() does not support dask with cupy backed DataArray.'), # noqa
)
out = mapper(agg)(agg.data, excludes)
return out
[docs]def mean(agg, passes=1, excludes=[np.nan], name='mean'):
"""
Returns Mean filtered array using a 3x3 window.
Default behaviour to 'mean' is to exclude NaNs from calculations.
Parameters
----------
agg : xarray.DataArray
2D array of input values to be filtered.
passes : int, default=1
Number of times to run mean.
name : str, default='mean'
Output xr.DataArray.name property.
Returns
-------
mean_agg : xarray.DataArray of same type as `agg`
2D aggregate array of filtered values.
Examples
--------
Focal mean works with NumPy backed xarray DataArray
.. sourcecode:: python
>>> import numpy as np
>>> import xarray as xr
>>> from xrspatial.focal import mean
>>> data = np.array([
[0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0.],
[0., 0., 9., 0., 0.],
[0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0.]])
>>> raster = xr.DataArray(data)
>>> mean_agg = mean(raster)
>>> print(mean_agg)
<xarray.DataArray 'mean' (dim_0: 5, dim_1: 5)>
array([[0., 0., 0., 0., 0.],
[0., 1., 1., 1., 0.],
[0., 1., 1., 1., 0.],
[0., 1., 1., 1., 0.],
[0., 0., 0., 0., 0.]])
Dimensions without coordinates: dim_0, dim_1
Focal mean works with Dask with NumPy backed xarray DataArray.
Increase number of runs by setting a specific value for parameter `passes`
.. sourcecode:: python
>>> import dask.array as da
>>> data_da = da.from_array(data, chunks=(3, 3))
>>> raster_da = xr.DataArray(data_da, dims=['y', 'x'], name='raster_da') # noqa
>>> print(raster_da)
<xarray.DataArray 'raster_da' (y: 5, x: 5)>
dask.array<array, shape=(5, 5), dtype=int64, chunksize=(3, 3), chunktype=numpy.ndarray> # noqa
Dimensions without coordinates: y, x
>>> mean_da = mean(raster_da, passes=2)
>>> print(mean_da)
<xarray.DataArray 'mean' (y: 5, x: 5)>
dask.array<_trim, shape=(5, 5), dtype=float64, chunksize=(3, 3), chunktype=numpy.ndarray> # noqa
Dimensions without coordinates: y, x
>>> print(mean_da.compute())
<xarray.DataArray 'mean' (y: 5, x: 5)>
array([[0.25 , 0.33333333, 0.5 , 0.33333333, 0.25 ],
[0.33333333, 0.44444444, 0.66666667, 0.44444444, 0.33333333],
[0.5 , 0.66666667, 1. , 0.66666667, 0.5 ],
[0.33333333, 0.44444444, 0.66666667, 0.44444444, 0.33333333],
[0.25 , 0.33333333, 0.5 , 0.33333333, 0.25 ]])
Dimensions without coordinates: y, x
Focal mean works with CuPy backed xarray DataArray.
In this example, we set `passes` to the number of elements of the array,
we'll get a mean array where every element has the same value.
.. sourcecode:: python
>>> import cupy
>>> raster_cupy = xr.DataArray(cupy.asarray(data), name='raster_cupy')
>>> mean_cupy = mean(raster_cupy, passes=25)
>>> print(type(mean_cupy.data))
<class 'cupy.core.core.ndarray'>
>>> print(mean_cupy)
<xarray.DataArray 'mean' (dim_0: 5, dim_1: 5)>
array([[0.47928995, 0.47928995, 0.47928995, 0.47928995, 0.47928995],
[0.47928995, 0.47928995, 0.47928995, 0.47928995, 0.47928995],
[0.47928995, 0.47928995, 0.47928995, 0.47928995, 0.47928995],
[0.47928995, 0.47928995, 0.47928995, 0.47928995, 0.47928995],
[0.47928995, 0.47928995, 0.47928995, 0.47928995, 0.47928995]])
Dimensions without coordinates: dim_0, dim_1
"""
out = agg.data.astype(float)
for i in range(passes):
out = _mean(out, tuple(excludes))
return DataArray(out,
name=name,
dims=agg.dims,
coords=agg.coords,
attrs=agg.attrs)
@ngjit
def _calc_mean(array):
return np.nanmean(array)
@ngjit
def _calc_sum(array):
return np.nansum(array)
@ngjit
def _calc_min(array):
return np.nanmin(array)
@ngjit
def _calc_max(array):
return np.nanmax(array)
@ngjit
def _calc_std(array):
return np.nanstd(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
def _apply_numpy(data, kernel, func):
data = data.astype(np.float32)
out = np.zeros_like(data)
rows, cols = data.shape
krows, kcols = kernel.shape
hrows, hcols = int(krows / 2), int(kcols / 2)
kernel_values = np.zeros_like(kernel, dtype=data.dtype)
for y in prange(rows):
for x in prange(cols):
# kernel values are all nans at the beginning of each step
kernel_values.fill(np.nan)
for ky in range(y - hrows, y + hrows + 1):
for kx in range(x - hcols, x + hcols + 1):
if ky >= 0 and ky < rows and kx >= 0 and kx < cols:
kyidx, kxidx = ky - (y - hrows), kx - (x - hcols)
if kernel[kyidx, kxidx] == 1:
kernel_values[kyidx, kxidx] = data[ky, kx]
out[y, x] = func(kernel_values)
return out
def _apply_dask_numpy(data, kernel, func):
data = data.astype(np.float32)
_func = partial(_apply_numpy, kernel=kernel, func=func)
pad_h = kernel.shape[0] // 2
pad_w = kernel.shape[1] // 2
out = data.map_overlap(_func,
depth=(pad_h, pad_w),
boundary=np.nan,
meta=np.array(()))
return out
[docs]def apply(raster, kernel, func=_calc_mean, name='focal_apply'):
"""
Returns custom function applied array using a user-created window.
Parameters
----------
raster : xarray.DataArray
2D array of input values to be filtered. Can be a NumPy backed,
or Dask with NumPy backed DataArray.
kernel : numpy.ndarray
2D array where values of 1 indicate the kernel.
func : callable, default=xrspatial.focal._calc_mean
Function which takes an input array and returns an array.
Returns
-------
agg : xarray.DataArray of same type as `raster`
2D aggregate array of filtered values.
Examples
--------
Focal apply works with NumPy backed xarray DataArray
.. sourcecode:: python
>>> import numpy as np
>>> import xarray as xr
>>> from xrspatial.convolution import circle_kernel
>>> from xrspatial.focal import apply
>>> data = np.arange(20, dtype=np.float64).reshape(4, 5)
>>> raster = xr.DataArray(data, dims=['y', 'x'], name='raster')
>>> print(raster)
<xarray.DataArray 'raster' (y: 4, x: 5)>
array([[ 0., 1., 2., 3., 4.],
[ 5., 6., 7., 8., 9.],
[10., 11., 12., 13., 14.],
[15., 16., 17., 18., 19.]])
Dimensions without coordinates: y, x
>>> kernel = circle_kernel(2, 2, 3)
>>> kernel
array([[0., 1., 0.],
[1., 1., 1.],
[0., 1., 0.]])
>>> # apply kernel mean by default
>>> apply_mean_agg = apply(raster, kernel)
>>> apply_mean_agg
<xarray.DataArray 'focal_apply' (y: 4, x: 5)>
array([[ 2. , 2.25 , 3.25 , 4.25 , 5.33333333],
[ 5.25 , 6. , 7. , 8. , 8.75 ],
[10.25 , 11. , 12. , 13. , 13.75 ],
[13.66666667, 14.75 , 15.75 , 16.75 , 17. ]])
Dimensions without coordinates: y, x
Focal apply works with Dask with NumPy backed xarray DataArray.
Note that if input raster is a numpy or dask with numpy backed data array,
the applied function must be decorated with ``numba.jit``
xrspatial already provides ``ngjit`` decorator, where:
``ngjit = numba.jit(nopython=True, nogil=True)``
.. sourcecode:: python
>>> from xrspatial.utils import ngjit
>>> from xrspatial.convolution import custom_kernel
>>> kernel = custom_kernel(np.array([
[0, 1, 0],
[0, 1, 1],
[0, 1, 0],
]))
>>> weight = np.array([
[0, 0.5, 0],
[0, 1, 0.5],
[0, 0.5, 0],
])
>>> @ngjit
>>> def func(kernel_data):
... weight = np.array([
... [0, 0.5, 0],
... [0, 1, 0.5],
... [0, 0.5, 0],
... ])
... return np.nansum(kernel_data * weight)
>>> import dask.array as da
>>> data_da = da.from_array(np.ones((6, 4), dtype=np.float64), chunks=(3, 2))
>>> raster_da = xr.DataArray(data_da, dims=['y', 'x'], name='raster_da')
>>> print(raster_da)
<xarray.DataArray 'raster_da' (y: 6, x: 4)>
dask.array<array, shape=(6, 4), dtype=float64, chunksize=(3, 2), chunktype=numpy.ndarray> # noqa
Dimensions without coordinates: y, x
>>> apply_func_agg = apply(raster_da, kernel, func)
>>> print(apply_func_agg)
<xarray.DataArray 'focal_apply' (y: 6, x: 4)>
dask.array<_trim, shape=(6, 4), dtype=float64, chunksize=(3, 2), chunktype=numpy.ndarray> # noqa
Dimensions without coordinates: y, x
>>> print(apply_func_agg.compute())
<xarray.DataArray 'focal_apply' (y: 6, x: 4)>
array([[2. , 2. , 2. , 1.5],
[2.5, 2.5, 2.5, 2. ],
[2.5, 2.5, 2.5, 2. ],
[2.5, 2.5, 2.5, 2. ],
[2.5, 2.5, 2.5, 2. ],
[2. , 2. , 2. , 1.5]])
Dimensions without coordinates: y, x
"""
# validate raster
if not isinstance(raster, DataArray):
raise TypeError("`raster` must be instance of DataArray")
if raster.ndim != 2:
raise ValueError("`raster` must be 2D")
# Validate the kernel
kernel = custom_kernel(kernel)
# apply kernel to raster values
# if raster is a numpy or dask with numpy backed data array,
# the function func must be a @ngjit
mapper = ArrayTypeFunctionMapping(
numpy_func=_apply_numpy,
cupy_func=lambda *args: not_implemented_func(
*args, messages='apply() does not support cupy backed DataArray.'),
dask_func=_apply_dask_numpy,
dask_cupy_func=lambda *args: not_implemented_func(
*args, messages='apply() does not support dask with cupy backed DataArray.'),
)
out = mapper(raster)(raster.data, kernel, func)
result = DataArray(out,
name=name,
coords=raster.coords,
dims=raster.dims,
attrs=raster.attrs)
return result
@cuda.jit
def _focal_min_cuda(data, kernel, out):
i, j = cuda.grid(2)
delta_rows = kernel.shape[0] // 2
delta_cols = kernel.shape[1] // 2
data_rows, data_cols = data.shape
if i < delta_rows or i >= data_rows - delta_rows or \
j < delta_cols or j >= data_cols - delta_cols:
return
s = data[i, j]
for k in range(kernel.shape[0]):
for h in range(kernel.shape[1]):
i_k = i + k - delta_rows
j_h = j + h - delta_cols
if (i_k >= 0) and (i_k < data_rows) and (j_h >= 0) and (j_h < data_cols):
if (kernel[k, h] != 0) and s > data[i_k, j_h]:
s = data[i_k, j_h]
out[i, j] = s
@cuda.jit
def _focal_max_cuda(data, kernel, out):
i, j = cuda.grid(2)
delta_rows = kernel.shape[0] // 2
delta_cols = kernel.shape[1] // 2
data_rows, data_cols = data.shape
if i < delta_rows or i >= data_rows - delta_rows or \
j < delta_cols or j >= data_cols - delta_cols:
return
s = data[i, j]
for k in range(kernel.shape[0]):
for h in range(kernel.shape[1]):
i_k = i + k - delta_rows
j_h = j + h - delta_cols
if (i_k >= 0) and (i_k < data_rows) and (j_h >= 0) and (j_h < data_cols):
if (kernel[k, h] != 0) and s < data[i_k, j_h]:
s = data[i_k, j_h]
out[i, j] = s
def _focal_range_cupy(data, kernel):
focal_min = _focal_stats_func_cupy(data, kernel, _focal_min_cuda)
focal_max = _focal_stats_func_cupy(data, kernel, _focal_max_cuda)
out = focal_max - focal_min
return out
@cuda.jit
def _focal_std_cuda(data, kernel, out):
i, j = cuda.grid(2)
delta_rows = kernel.shape[0] // 2
delta_cols = kernel.shape[1] // 2
data_rows, data_cols = data.shape
if i < delta_rows or i >= data_rows - delta_rows or \
j < delta_cols or j >= data_cols - delta_cols:
return
sum_squares = 0
sum = 0
count = 0
for k in range(kernel.shape[0]):
for h in range(kernel.shape[1]):
i_k = i + k - delta_rows
j_h = j + h - delta_cols
if (i_k >= 0) and (i_k < data_rows) and (j_h >= 0) and (j_h < data_cols):
sum_squares += (kernel[k, h]*data[i_k, j_h])**2
sum += kernel[k, h]*data[i_k, j_h]
count += kernel[k, h]
squared_sum = sum**2
out[i, j] = sqrt((sum_squares - squared_sum/count) / count)
@cuda.jit
def _focal_var_cuda(data, kernel, out):
i, j = cuda.grid(2)
delta_rows = kernel.shape[0] // 2
delta_cols = kernel.shape[1] // 2
data_rows, data_cols = data.shape
if i < delta_rows or i >= data_rows - delta_rows or \
j < delta_cols or j >= data_cols - delta_cols:
return
sum_squares = 0
sum = 0
count = 0
for k in range(kernel.shape[0]):
for h in range(kernel.shape[1]):
i_k = i + k - delta_rows
j_h = j + h - delta_cols
if (i_k >= 0) and (i_k < data_rows) and (j_h >= 0) and (j_h < data_cols):
sum_squares += (kernel[k, h]*data[i_k, j_h])**2
sum += kernel[k, h]*data[i_k, j_h]
count += kernel[k, h]
squared_sum = sum**2
out[i, j] = (sum_squares - squared_sum/count) / count
def _focal_mean_cupy(data, kernel):
out = convolve_2d(data, kernel / kernel.sum())
return out
def _focal_sum_cupy(data, kernel):
out = convolve_2d(data, kernel)
return out
def _focal_stats_func_cupy(data, kernel, func=_focal_max_cuda):
out = cupy.empty(data.shape, dtype='f4')
out[:, :] = cupy.nan
griddim, blockdim = cuda_args(data.shape)
func[griddim, blockdim](data, kernel, cupy.asarray(out))
return out
def _focal_stats_cupy(agg, kernel, stats_funcs):
_stats_cupy_mapper = dict(
mean=_focal_mean_cupy,
sum=_focal_sum_cupy,
range=_focal_range_cupy,
max=lambda *args: _focal_stats_func_cupy(*args, func=_focal_max_cuda),
min=lambda *args: _focal_stats_func_cupy(*args, func=_focal_min_cuda),
std=lambda *args: _focal_stats_func_cupy(*args, func=_focal_std_cuda),
var=lambda *args: _focal_stats_func_cupy(*args, func=_focal_var_cuda),
)
stats_aggs = []
for stats in stats_funcs:
data = agg.data.astype(cupy.float32)
stats_data = _stats_cupy_mapper[stats](data, kernel)
stats_agg = xr.DataArray(
stats_data,
dims=agg.dims,
coords=agg.coords,
attrs=agg.attrs
)
stats_aggs.append(stats_agg)
stats = xr.concat(stats_aggs, pd.Index(stats_funcs, name='stats'))
return stats
def _focal_stats_cpu(agg, kernel, stats_funcs):
_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
[docs]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. Can be a NumPy backed,
Cupy backed, or Dask with NumPy backed DataArray.
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.
Examples
--------
.. sourcecode:: python
>>> import numpy as np
>>> import xarray as xr
>>> from xrspatial.convolution import circle_kernel
>>> kernel = circle_kernel(1, 1, 1)
>>> kernel
array([[0., 1., 0.],
[1., 1., 1.],
[0., 1., 0.]])
>>> data = np.array([
[0, 0, 0, 0, 0, 0],
[1, 1, 2, 2, 1, 1],
[2, 2, 1, 1, 2, 2],
[3, 3, 0, 0, 3, 3],
])
>>> from xrspatial.focal import focal_stats
>>> focal_stats(xr.DataArray(data), kernel, stats_funcs=['min', 'sum'])
<xarray.DataArray 'focal_apply' (stats: 2, dim_0: 4, dim_1: 6)>
array([[[0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0.],
[1., 1., 0., 0., 1., 1.],
[2., 0., 0., 0., 0., 2.]],
[[1., 1., 2., 2., 1., 1.],
[4., 6., 6., 6., 6., 4.],
[8., 9., 6., 6., 9., 8.],
[8., 8., 4., 4., 8., 8.]]])
Coordinates:
* stats (stats) object 'min' 'sum'
Dimensions without coordinates: dim_0, dim_1
"""
# validate raster
if not isinstance(agg, DataArray):
raise TypeError("`agg` must be instance of DataArray")
if agg.ndim != 2:
raise ValueError("`agg` must be 2D")
# Validate the kernel
kernel = custom_kernel(kernel)
mapper = ArrayTypeFunctionMapping(
numpy_func=_focal_stats_cpu,
cupy_func=_focal_stats_cupy,
dask_func=_focal_stats_cpu,
dask_cupy_func=lambda *args: not_implemented_func(
*args, messages='focal_stats() does not support dask with cupy backed DataArray.'),
)
result = mapper(agg)(agg, kernel, stats_funcs)
return result
@ngjit
def _calc_hotspots_numpy(z_array):
out = np.zeros_like(z_array, dtype=np.int8)
rows, cols = z_array.shape
for y in prange(rows):
for x in prange(cols):
zscore = z_array[y, x]
# find p value
p_value = 1.0
if abs(zscore) >= 2.33:
p_value = 0.0099
elif abs(zscore) >= 1.65:
p_value = 0.0495
elif abs(zscore) >= 1.29:
p_value = 0.0985
# confidence
confidence = 0
if abs(zscore) > 2.58 and p_value < 0.01:
confidence = 99
elif abs(zscore) > 1.96 and p_value < 0.05:
confidence = 95
elif abs(zscore) > 1.65 and p_value < 0.1:
confidence = 90
hot_cold = 0
if zscore > 0:
hot_cold = 1
elif zscore < 0:
hot_cold = -1
out[y, x] = hot_cold * confidence
return out
def _hotspots_numpy(raster, kernel):
if not (issubclass(raster.data.dtype.type, np.integer) or
issubclass(raster.data.dtype.type, np.floating)):
raise ValueError("data type must be integer or float")
data = raster.data.astype(np.float32)
# apply kernel to raster values
mean_array = convolve_2d(data, kernel / kernel.sum())
# calculate z-scores
global_mean = np.nanmean(data)
global_std = np.nanstd(data)
if global_std == 0:
raise ZeroDivisionError(
"Standard deviation of the input raster values is 0."
)
z_array = (mean_array - global_mean) / global_std
out = _calc_hotspots_numpy(z_array)
return out
def _hotspots_dask_numpy(raster, kernel):
data = raster.data.astype(np.float32)
# apply kernel to raster values
mean_array = convolve_2d(data, kernel / kernel.sum())
# calculate z-scores
global_mean = da.nanmean(data)
global_std = da.nanstd(data)
# commented out to avoid early compute to check if global_std is zero
# if global_std == 0:
# raise ZeroDivisionError(
# "Standard deviation of the input raster values is 0."
# )
z_array = (mean_array - global_mean) / global_std
_func = partial(_calc_hotspots_numpy)
pad_h = kernel.shape[0] // 2
pad_w = kernel.shape[1] // 2
out = z_array.map_overlap(_func,
depth=(pad_h, pad_w),
boundary=np.nan,
meta=np.array(()))
return out
@nb.cuda.jit(device=True)
def _gpu_hotspots(data):
zscore = data[0, 0]
# find p value
p_value = 1.0
if abs(zscore) >= 2.33:
p_value = 0.0099
elif abs(zscore) >= 1.65:
p_value = 0.0495
elif abs(zscore) >= 1.29:
p_value = 0.0985
# confidence
confidence = 0
if abs(zscore) > 2.58 and p_value < 0.01:
confidence = 99
elif abs(zscore) > 1.96 and p_value < 0.05:
confidence = 95
elif abs(zscore) > 1.65 and p_value < 0.1:
confidence = 90
hot_cold = 0
if zscore > 0:
hot_cold = 1
elif zscore < 0:
hot_cold = -1
return hot_cold * confidence
@nb.cuda.jit
def _run_gpu_hotspots(data, out):
i, j = nb.cuda.grid(2)
if i >= 0 and i < out.shape[0] and j >= 0 and j < out.shape[1]:
out[i, j] = _gpu_hotspots(data[i:i + 1, j:j + 1])
def _hotspots_cupy(raster, kernel):
if not (issubclass(raster.data.dtype.type, cupy.integer) or
issubclass(raster.data.dtype.type, cupy.floating)):
raise ValueError("data type must be integer or float")
data = raster.data.astype(cupy.float32)
# apply kernel to raster values
mean_array = convolve_2d(data, kernel / kernel.sum())
# calculate z-scores
global_mean = cupy.nanmean(data)
global_std = cupy.nanstd(data)
if global_std == 0:
raise ZeroDivisionError(
"Standard deviation of the input raster values is 0."
)
z_array = (mean_array - global_mean) / global_std
out = cupy.zeros_like(z_array, dtype=cupy.int8)
griddim, blockdim = cuda_args(z_array.shape)
_run_gpu_hotspots[griddim, blockdim](z_array, out)
return out
[docs]def hotspots(raster, kernel):
"""
Identify statistically significant hot spots and cold spots in an
input raster. To be a statistically significant hot spot, a feature
will have a high value and be surrounded by other features with
high values as well.
Neighborhood of a feature defined by the input kernel, which
currently support a shape of circle, annulus, or custom kernel.
The result should be a raster with the following 7 values:
- 90 for 90% confidence high value cluster
- 95 for 95% confidence high value cluster
- 99 for 99% confidence high value cluster
- 90 for 90% confidence low value cluster
- 95 for 95% confidence low value cluster
- 99 for 99% confidence low value cluster
- 0 for no significance
Parameters
----------
raster : xarray.DataArray
2D Input raster image with `raster.shape` = (height, width).
Can be a NumPy backed, CuPy backed, or Dask with NumPy backed DataArray
kernel : Numpy Array
2D array where values of 1 indicate the kernel.
Returns
-------
hotspots_agg : xarray.DataArray of same type as `raster`
2D array of hotspots with values indicating confidence level.
Examples
--------
.. sourcecode:: python
>>> import numpy as np
>>> import xarray as xr
>>> from xrspatial.convolution import custom_kernel
>>> kernel = custom_kernel(np.array([[1, 1, 0]]))
>>> data = np.array([
... [0, 1000, 1000, 0, 0, 0],
... [0, 0, 0, -1000, -1000, 0],
... [0, -900, -900, 0, 0, 0],
... [0, 100, 1000, 0, 0, 0]])
>>> from xrspatial.focal import hotspots
>>> hotspots(xr.DataArray(data), kernel)
array([[ 0, 0, 95, 0, 0, 0],
[ 0, 0, 0, 0, -90, 0],
[ 0, 0, -90, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0]], dtype=int8)
Dimensions without coordinates: dim_0, dim_1
"""
# validate raster
if not isinstance(raster, DataArray):
raise TypeError("`raster` must be instance of DataArray")
if raster.ndim != 2:
raise ValueError("`raster` must be 2D")
mapper = ArrayTypeFunctionMapping(
numpy_func=_hotspots_numpy,
cupy_func=_hotspots_cupy,
dask_func=_hotspots_dask_numpy,
dask_cupy_func=lambda *args: not_implemented_func(
*args, messages='hotspots() does not support dask with cupy backed DataArray.'), # noqa
)
out = mapper(raster)(raster, kernel)
attrs = copy.deepcopy(raster.attrs)
attrs['unit'] = '%'
return DataArray(out,
coords=raster.coords,
dims=raster.dims,
attrs=attrs)