Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
317 changes: 141 additions & 176 deletions xrspatial/convolution.py

Large diffs are not rendered by default.

183 changes: 65 additions & 118 deletions xrspatial/focal.py
Original file line number Diff line number Diff line change
Expand Up @@ -173,46 +173,31 @@ def mean(agg, passes=1, excludes=[np.nan], name='mean'):
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' (dim_0: 5, dim_1: 5)>
<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: dim_0, dim_1
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
>>> print(mean(raster, passes=25))
>>> 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.47928994, 0.47928994, 0.47928994, 0.47928994, 0.47928994],
[0.47928994, 0.47928994, 0.47928994, 0.47928994, 0.47928994],
[0.47928994, 0.47928994, 0.47928994, 0.47928994, 0.47928994],
[0.47928994, 0.47928994, 0.47928994, 0.47928994, 0.47928994],
[0.47928994, 0.47928994, 0.47928994, 0.47928994, 0.47928994]])
Dimensions without coordinates: dim_0, dim_1

.. sourcecode:: python
>>> import cupy
>>> data_cupy = cupy.asarray([
[0, 1, 1, 1, 1, 2],
[0, 0, 1, 1, 2, 2],
[0, -1, 0, 2, 2, 2],
[-2, -2, -1, 0, 1, 1],
])
>>> raster_cupy = xr.DataArray(data_cupy, dims=['y', 'x'])
>>> mean_cupy = mean(raster_cupy)
>>> print(type(mean_cupy.data))
<class 'cupy.core.core.ndarray'>
>>> print(mean_cupy)
<xarray.DataArray 'mean' (y: 4, x: 6)>
array([[0.25 , 0.5 , 0.8333333 , 1.1666666, 1.5 , 1.75 ], # noqa
[0. , 0.22222222, 0.6666667 , 1.2222222, 1.6666666, 1.8333334], # noqa
[-0.8333333, -0.5555556, 0. , 0.8888889, 1.4444444, 1.6666666], # noqa
[-1.25 ,-1. , -0.33333334, 0.6666667, 1.3333334, 1.5 ]], dtype=float32) # noqa
Dimensions without coordinates: y, x
"""

out = agg.data
Expand Down Expand Up @@ -331,8 +316,9 @@ def apply(raster, kernel, func=_calc_mean, name='focal_apply'):
Parameters
----------
raster : xarray.DataArray
2D array of input values to be filtered.
kernel : numpy.array
2D array of input values to be filtered. Can be a NumPy backed,
CuPy 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.
Expand Down Expand Up @@ -463,7 +449,8 @@ def focal_stats(agg,
Parameters
----------
agg : xarray.DataArray
2D array of input values to be analysed.
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
Expand All @@ -475,6 +462,38 @@ def focal_stats(agg,
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
"""

_function_mapping = {
Expand Down Expand Up @@ -661,6 +680,7 @@ def hotspots(raster, kernel):
----------
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.

Expand All @@ -671,100 +691,27 @@ def hotspots(raster, kernel):

Examples
--------
.. plot::
:include-source:

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

from xrspatial import generate_terrain, aspect
from xrspatial.convolution import circle_kernel
from xrspatial.focal import hotspots


# Generate Example Terrain
W = 500
H = 300

template_terrain = xr.DataArray(np.zeros((H, W)))
x_range=(-20e6, 20e6)
y_range=(-20e6, 20e6)

terrain_agg = generate_terrain(
template_terrain, x_range=x_range, y_range=y_range
)

# Edit Attributes
terrain_agg = terrain_agg.assign_attrs(
{
'Description': 'Example Terrain',
'units': 'km',
'Max Elevation': '4000',
}
)

terrain_agg = terrain_agg.rename({'x': 'lon', 'y': 'lat'})
terrain_agg = terrain_agg.rename('Elevation')

# Create Kernel
kernel = circle_kernel(10, 10, 100)

# Create Hotspots Aggregate array
hotspots_agg = hotspots(raster = terrain_agg,
kernel = kernel)

# Edit Attributes
hotspots_agg = hotspots_agg.rename('Significance')
hotspots_agg = hotspots_agg.assign_attrs(
{
'Description': 'Example Hotspots',
'units': '%',
}
)

# Plot Terrain
terrain_agg.plot(cmap = 'terrain', aspect = 2, size = 4)
plt.title("Terrain")
plt.ylabel("latitude")
plt.xlabel("longitude")

# Plot Hotspots
hotspots_agg.plot(aspect = 2, size = 4)
plt.title("Hotspots")
plt.ylabel("latitude")
plt.xlabel("longitude")

.. sourcecode:: python

>>> print(terrain_agg[200:203, 200:202])
<xarray.DataArray 'Elevation' (lat: 3, lon: 2)>
array([[1264.02296597, 1261.947921 ],
[1285.37105519, 1282.48079719],
[1306.02339636, 1303.4069579 ]])
Coordinates:
* lon (lon) float64 -3.96e+06 -3.88e+06
* lat (lat) float64 6.733e+06 6.867e+06 7e+06
Attributes:
res: (80000.0, 133333.3333333333)
Description: Example Terrain
units: km
Max Elevation: 4000

>>> print(hotspots_agg[200:203, 200:202])
<xarray.DataArray 'Significance' (lat: 3, lon: 2)>
array([[0, 0],
[0, 0],
[0, 0]], dtype=int8)
Coordinates:
* lon (lon) float64 -3.96e+06 -3.88e+06
* lat (lat) float64 6.733e+06 6.867e+06 7e+06
Attributes:
res: (80000.0, 133333.3333333333)
Description: Example Hotspots
units: %
Max Elevation: 4000
>>> 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
"""

# TODO: edit unit of output raster to percent (%)

# validate raster
if not isinstance(raster, DataArray):
raise TypeError("`raster` must be instance of DataArray")
Expand Down
36 changes: 14 additions & 22 deletions xrspatial/hillshade.py
Original file line number Diff line number Diff line change
Expand Up @@ -128,29 +128,21 @@ def hillshade(agg: xr.DataArray,

Examples
--------
.. plot::
:include-source:

import numpy as np
import xarray as xr
from xrspatial import hillshade

data = np.array([
[0., 0., 0., 0., 0.],
[0., 1., 0., 2., 0.],
[0., 0., 3., 0., 0.],
[0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0.]
])
n, m = data.shape
raster = xr.DataArray(data, dims=['y', 'x'], name='raster')
raster['y'] = np.arange(n)[::-1]
raster['x'] = np.arange(m)

hillshade_agg = hillshade(raster)

.. sourcecode:: python

>>> import numpy as np
>>> import xarray as xr
>>> from xrspatial import hillshade
>>> data = np.array([
... [0., 0., 0., 0., 0.],
... [0., 1., 0., 2., 0.],
... [0., 0., 3., 0., 0.],
... [0., 0., 0., 0., 0.],
... [0., 0., 0., 0., 0.]])
>>> n, m = data.shape
>>> raster = xr.DataArray(data, dims=['y', 'x'], name='raster')
>>> raster['y'] = np.arange(n)[::-1]
>>> raster['x'] = np.arange(m)
>>> hillshade_agg = hillshade(raster)
>>> print(hillshade_agg)
<xarray.DataArray 'hillshade' (y: 5, x: 5)>
array([[ nan, nan, nan, nan, nan],
Expand Down
Loading