Skip to content

Fixes #906: fill remaining dask+cupy gaps (terrain, perlin, crosstab) #913

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
brendancol merged 4 commits into from
Feb 27, 2026
Merged

Fixes #906: fill remaining dask+cupy gaps (terrain, perlin, crosstab) #913

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5f6c332
Fixes #906: add dask+cupy backends for perlin, terrain, crosstab
brendancol Feb 27, 2026
314a65e
Update README feature matrix for new and existing GPU backends
brendancol Feb 27, 2026
6c6fa8a
Add dask/cupy/dask+cupy support for trim() and crop()
brendancol Feb 27, 2026
1455f45
Add native GPU kernel for zonal.apply() with automatic CPU fallback
brendancol Feb 27, 2026
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
14 changes: 7 additions & 7 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -224,10 +224,10 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
| [Curvature](xrspatial/curvature.py) | Measures rate of slope change (concavity/convexity) at each cell | ✅️ |✅️ |✅️ | ✅️ |
| [Hillshade](xrspatial/hillshade.py) | Simulates terrain illumination from a given sun angle and azimuth | ✅️ | ✅️ | ✅️ | ✅️ |
| [Slope](xrspatial/slope.py) | Computes terrain gradient steepness at each cell in degrees | ✅️ | ✅️ | ✅️ | ✅️ |
| [Terrain Generation](xrspatial/terrain.py) | Generates synthetic terrain elevation using fractal noise | ✅️ | ✅️ | ✅️ | |
| [Terrain Generation](xrspatial/terrain.py) | Generates synthetic terrain elevation using fractal noise | ✅️ | ✅️ | ✅️ | ✅️ |
| [Viewshed](xrspatial/viewshed.py) | Determines visible cells from a given observer point on terrain | ✅️ | ✅️ | ✅️ | ✅️ |
| [Min Observable Height](xrspatial/experimental/min_observable_height.py) | Finds the minimum observer height needed to see each cell *(experimental)* | ✅️ | | | |
| [Perlin Noise](xrspatial/perlin.py) | Generates smooth continuous random noise for procedural textures | ✅️ | ✅️ | ✅️ | |
| [Perlin Noise](xrspatial/perlin.py) | Generates smooth continuous random noise for procedural textures | ✅️ | ✅️ | ✅️ | ✅️ |
| [Bump Mapping](xrspatial/bump.py) | Adds randomized bump features to simulate natural terrain variation | ✅️ | | | |

-----------
Expand All @@ -236,12 +236,12 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e

| Name | Description | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
|:----------:|:------------|:----------------------:|:--------------------:|:-------------------:|:------:|
| [Apply](xrspatial/zonal.py) | Applies a custom function to each zone in a classified raster | ✅️ | ✅️ | | |
| [Crop](xrspatial/zonal.py) | Extracts the bounding rectangle of a specific zone | ✅️ | | | |
| [Apply](xrspatial/zonal.py) | Applies a custom function to each zone in a classified raster | ✅️ | ✅️ | ✅️ | ✅️ |
| [Crop](xrspatial/zonal.py) | Extracts the bounding rectangle of a specific zone | ✅️ | ✅️ | ✅️ | ✅️ |
| [Regions](xrspatial/zonal.py) | Identifies connected regions of non-zero cells | ✅️ | ✅️ | ✅️ | ✅️ |
| [Trim](xrspatial/zonal.py) | Removes nodata border rows and columns from a raster | ✅️ | | | |
| [Zonal Statistics](xrspatial/zonal.py) | Computes summary statistics for a value raster within each zone | ✅️ | ✅️| | |
| [Zonal Cross Tabulate](xrspatial/zonal.py) | Cross-tabulates agreement between two categorical rasters | ✅️ | ✅️| | |
| [Trim](xrspatial/zonal.py) | Removes nodata border rows and columns from a raster | ✅️ | ✅️ | ✅️ | ✅️ |
| [Zonal Statistics](xrspatial/zonal.py) | Computes summary statistics for a value raster within each zone | ✅️ | ✅️| ✅️ | 🔄 |
| [Zonal Cross Tabulate](xrspatial/zonal.py) | Cross-tabulates agreement between two categorical rasters | ✅️ | ✅️| 🔄 | 🔄 |

#### Usage

Expand Down
35 changes: 32 additions & 3 deletions xrspatial/perlin.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -192,6 +192,37 @@ def _perlin_cupy(data: cupy.ndarray,
return data


def _perlin_dask_cupy(data: da.Array,
freq: tuple,
seed: int) -> da.Array:
np.random.seed(seed)
p = cupy.asarray(np.random.permutation(2**20))
p = cupy.append(p, p)

height, width = data.shape

def _chunk_perlin(block, block_info=None):
info = block_info[0]
y_start, y_end = info['array-location'][0]
x_start, x_end = info['array-location'][1]
x0 = freq[0] * x_start / width
x1 = freq[0] * x_end / width
y0 = freq[1] * y_start / height
y1 = freq[1] * y_end / height

out = cupy.empty(block.shape, dtype=cupy.float32)
griddim, blockdim = cuda_args(block.shape)
_perlin_gpu[griddim, blockdim](p, x0, x1, y0, y1, 1.0, out)
return out

data = da.map_blocks(_chunk_perlin, data, dtype=cupy.float32,
meta=cupy.array((), dtype=cupy.float32))

min_val, max_val = dask.compute(da.min(data), da.max(data))
data = (data - min_val) / (max_val - min_val)
return data


def perlin(agg: xr.DataArray,
freq: tuple = (1, 1),
seed: int = 5,
Expand Down Expand Up @@ -246,9 +277,7 @@ def perlin(agg: xr.DataArray,
numpy_func=_perlin_numpy,
cupy_func=_perlin_cupy,
dask_func=_perlin_dask_numpy,
dask_cupy_func=lambda *args: not_implemented_func(
*args, messages='perlin() does not support dask with cupy backed DataArray', # noqa
)
dask_cupy_func=_perlin_dask_cupy
)
out = mapper(agg)(agg.data, freq, seed)
result = xr.DataArray(out,
Expand Down
53 changes: 50 additions & 3 deletions xrspatial/terrain.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,55 @@ def _terrain_gpu(height_map, seed, x_range=(0, 1), y_range=(0, 1)):
return height_map


def _terrain_dask_cupy(data: da.Array,
seed: int,
x_range_scaled: tuple,
y_range_scaled: tuple,
zfactor: int) -> da.Array:
data = data * 0
height, width = data.shape

NOISE_LAYERS = ((1 / 2**i, (2**i, 2**i)) for i in range(16))

nrange = np.arange(2**20, dtype=np.int32)

for i, (m, (xfreq, yfreq)) in enumerate(NOISE_LAYERS):
np.random.seed(seed + i)
p = cupy.asarray(np.random.permutation(nrange))
p = cupy.append(p, p)

def _chunk_noise(block, block_info=None, _p=p, _m=m,
_xr=x_range_scaled, _yr=y_range_scaled,
_xf=xfreq, _yf=yfreq):
info = block_info[0]
y_start, y_end = info['array-location'][0]
x_start, x_end = info['array-location'][1]
x0 = _xr[0] + (_xr[1] - _xr[0]) * x_start / width
x1 = _xr[0] + (_xr[1] - _xr[0]) * x_end / width
y0 = _yr[0] + (_yr[1] - _yr[0]) * y_start / height
y1 = _yr[0] + (_yr[1] - _yr[0]) * y_end / height

out = cupy.empty(block.shape, dtype=cupy.float32)
griddim, blockdim = cuda_args(block.shape)
_perlin_gpu[griddim, blockdim](
_p, x0 * _xf, x1 * _xf, y0 * _yf, y1 * _yf, _m, out
)
return out

noise = da.map_blocks(_chunk_noise, data, dtype=cupy.float32,
meta=cupy.array((), dtype=cupy.float32))
data = data + noise

data /= (1.00 + 0.50 + 0.25 + 0.13 + 0.06 + 0.03)
data = data ** 3

min_val, max_val = dask.compute(da.min(data), da.max(data))
data = (data - min_val) / (max_val - min_val)
data = da.where(data < 0.3, 0, data)
data *= zfactor
return data


def _terrain_cupy(data: cupy.ndarray,
seed: int,
x_range_scaled: tuple,
Expand Down Expand Up @@ -263,9 +312,7 @@ def generate_terrain(agg: xr.DataArray,
numpy_func=_terrain_numpy,
cupy_func=_terrain_cupy,
dask_func=_terrain_dask_numpy,
dask_cupy_func=lambda *args: not_implemented_func(
*args, messages='generate_terrain() does not support dask with cupy backed DataArray' # noqa
)
dask_cupy_func=_terrain_dask_cupy
)
out = mapper(agg)(agg.data, seed, x_range_scaled, y_range_scaled, zfactor)
canvas = ds.Canvas(
Expand Down
Loading