Fixes #903: add comprehensive input validation across public API

Add _validate_raster() and _validate_scalar() utilities to utils.py and
call them at the top of every public function (~35 entry points across
15 modules). Invalid input now raises clear TypeError/ValueError with
function name, parameter name, and expected-vs-actual values instead of
cryptic numba TypingError deep in JIT kernels.

- _validate_raster: checks DataArray type, ndim, numeric/integer dtype
- _validate_scalar: checks scalar type (incl. numpy equivalents), range
- Replace ad-hoc isinstance/ndim checks in focal, zonal, cost_distance
- Add validation to all surface, proximity, classify, multispectral,
  viewshed, perlin, terrain, and bump functions
- Add test_validation.py with 74 systematic tests covering type, ndim,
  dtype, scalar range, and integer-only checks

Author: brendancol

xrspatial/aspect.py

@@ -20,6 +20,7 @@
 from xrspatial.utils import _extract_latlon_coords
 from xrspatial.utils import _pad_array
 from xrspatial.utils import _validate_boundary
+from xrspatial.utils import _validate_raster
 from xrspatial.utils import cuda_args
 from xrspatial.utils import ngjit
 from xrspatial.dataset_support import supports_dataset
@@ -396,6 +397,8 @@ def aspect(agg: xr.DataArray,
         >>> aspect_agg = aspect(raster)
     """
 
+    _validate_raster(agg, func_name='aspect', name='agg')
+
     if method not in ('planar', 'geodesic'):
         raise ValueError(
             f"method must be 'planar' or 'geodesic', got {method!r}"

xrspatial/bump.py

@@ -4,7 +4,7 @@
 import xarray as xr
 from xarray import DataArray
 
-from xrspatial.utils import ngjit
+from xrspatial.utils import _validate_scalar, ngjit
 
 # TODO: change parameters to take agg instead of height / width
 
@@ -194,6 +194,9 @@ def heights(locations, src, src_range, height = 20):
             Description:  Example Bump Map
             units:        km
     """
+    _validate_scalar(width, func_name='bump', name='width', dtype=int, min_val=1)
+    _validate_scalar(height, func_name='bump', name='height', dtype=int, min_val=1)
+
     linx = range(width)
     liny = range(height)
 

xrspatial/classify.py

@@ -24,7 +24,14 @@ class cupy(object):
 import numba as nb
 import numpy as np
 
-from xrspatial.utils import ArrayTypeFunctionMapping, cuda_args, ngjit, not_implemented_func
+from xrspatial.utils import (
+    ArrayTypeFunctionMapping,
+    _validate_raster,
+    _validate_scalar,
+    cuda_args,
+    ngjit,
+    not_implemented_func,
+)
 from xrspatial.dataset_support import supports_dataset
 
 
@@ -136,6 +143,7 @@ def binary(agg, values, name='binary'):
                [0.,  0.,  0.,  0.,  np.nan]], dtype=float32)
         Dimensions without coordinates: dim_0, dim_1
     """
+    _validate_raster(agg, func_name='binary', name='agg', ndim=None)
 
     mapper = ArrayTypeFunctionMapping(numpy_func=_run_numpy_binary,
                                       dask_func=_run_dask_numpy_binary,
@@ -380,6 +388,7 @@ def reclassify(agg: xr.DataArray,
 
     Reclassify works with Dask with CuPy backed xarray DataArray.
     """
+    _validate_raster(agg, func_name='reclassify', name='agg', ndim=None)
 
     if len(bins) != len(new_values):
         raise ValueError(
@@ -515,6 +524,8 @@ def quantile(agg: xr.DataArray,
         Attributes:
             res:      (10.0, 10.0)
     """
+    _validate_raster(agg, func_name='quantile', name='agg', ndim=None)
+    _validate_scalar(k, func_name='quantile', name='k', dtype=int, min_val=2)
 
     q = _quantile(agg, num_sample, k)
     k_q = q.shape[0]
@@ -836,6 +847,8 @@ def natural_breaks(agg: xr.DataArray,
                [ 4.,  4.,  4.,  4., nan]], dtype=float32)
         Dimensions without coordinates: dim_0, dim_1
     """
+    _validate_raster(agg, func_name='natural_breaks', name='agg', ndim=None)
+    _validate_scalar(k, func_name='natural_breaks', name='k', dtype=int, min_val=2)
 
     mapper = ArrayTypeFunctionMapping(
         numpy_func=_run_natural_break,
@@ -944,6 +957,8 @@ def equal_interval(agg: xr.DataArray,
         Attributes:
             res:      (10.0, 10.0)
     """
+    _validate_raster(agg, func_name='equal_interval', name='agg', ndim=None)
+    _validate_scalar(k, func_name='equal_interval', name='k', dtype=int, min_val=1)
 
     mapper = ArrayTypeFunctionMapping(
         numpy_func=lambda *args: _run_equal_interval(*args, module=np),
@@ -1015,6 +1030,8 @@ def std_mean(agg: xr.DataArray,
     ----------
         - PySAL: https://pysal.org/mapclassify/_modules/mapclassify/classifiers.html#StdMean
     """
+    _validate_raster(agg, func_name='std_mean', name='agg', ndim=None)
+
     mapper = ArrayTypeFunctionMapping(
         numpy_func=lambda *args: _run_std_mean(*args, module=np),
         dask_func=lambda *args: _run_std_mean(*args, module=da),
@@ -1112,6 +1129,8 @@ def head_tail_breaks(agg: xr.DataArray,
     ----------
         - PySAL: https://pysal.org/mapclassify/_modules/mapclassify/classifiers.html#HeadTailBreaks
     """
+    _validate_raster(agg, func_name='head_tail_breaks', name='agg', ndim=None)
+
     mapper = ArrayTypeFunctionMapping(
         numpy_func=lambda *args: _run_head_tail_breaks(*args, module=np),
         dask_func=_run_dask_head_tail_breaks,
@@ -1191,6 +1210,8 @@ def percentiles(agg: xr.DataArray,
     ----------
         - PySAL: https://pysal.org/mapclassify/_modules/mapclassify/classifiers.html#Percentiles
     """
+    _validate_raster(agg, func_name='percentiles', name='agg', ndim=None)
+
     if pct is None:
         pct = [1, 10, 50, 90, 99]
 
@@ -1328,6 +1349,9 @@ def maximum_breaks(agg: xr.DataArray,
     ----------
         - PySAL: https://pysal.org/mapclassify/_modules/mapclassify/classifiers.html#MaximumBreaks
     """
+    _validate_raster(agg, func_name='maximum_breaks', name='agg', ndim=None)
+    _validate_scalar(k, func_name='maximum_breaks', name='k', dtype=int, min_val=2)
+
     mapper = ArrayTypeFunctionMapping(
         numpy_func=lambda *args: _run_maximum_breaks(*args, module=np),
         dask_func=_run_dask_maximum_breaks,
@@ -1431,6 +1455,8 @@ def box_plot(agg: xr.DataArray,
     ----------
         - PySAL: https://pysal.org/mapclassify/_modules/mapclassify/classifiers.html#BoxPlot
     """
+    _validate_raster(agg, func_name='box_plot', name='agg', ndim=None)
+
     mapper = ArrayTypeFunctionMapping(
         numpy_func=lambda *args: _run_box_plot(*args, module=np),
         dask_func=lambda *args: _run_box_plot(*args, module=da),

xrspatial/cost_distance.py

@@ -49,6 +49,7 @@ class cupy:  # type: ignore[no-redef]
         ndarray = False
 
 from xrspatial.utils import (
+    _validate_raster,
     cuda_args, get_dataarray_resolution, ngjit,
     has_cuda_and_cupy, is_cupy_array, is_dask_cupy,
 )
@@ -1098,10 +1099,8 @@ def cost_distance(
         Source pixels have cost 0.  Unreachable pixels are NaN.
     """
     # --- validation ---
-    if raster.ndim != 2:
-        raise ValueError("raster must be 2-D")
-    if friction.ndim != 2:
-        raise ValueError("friction must be 2-D")
+    _validate_raster(raster, func_name='cost_distance', name='raster')
+    _validate_raster(friction, func_name='cost_distance', name='friction')
     if raster.shape != friction.shape:
         raise ValueError("raster and friction must have the same shape")
     if raster.dims != (y, x):

xrspatial/curvature.py

@@ -25,6 +25,7 @@ class cupy(object):
 from xrspatial.utils import _boundary_to_dask
 from xrspatial.utils import _pad_array
 from xrspatial.utils import _validate_boundary
+from xrspatial.utils import _validate_raster
 from xrspatial.utils import cuda_args
 from xrspatial.utils import get_dataarray_resolution
 from xrspatial.utils import ngjit
@@ -249,6 +250,8 @@ def curvature(agg: xr.DataArray,
         Attributes:
             res:      (10, 10)
     """
+    _validate_raster(agg, func_name='curvature', name='agg')
+
     cellsize_x, cellsize_y = get_dataarray_resolution(agg)
     cellsize = (cellsize_x + cellsize_y) / 2
 

xrspatial/focal.py

@@ -29,7 +29,8 @@ class cupy(object):
 
 from xrspatial.convolution import convolve_2d, custom_kernel, _convolve_2d_numpy
 from xrspatial.utils import (ArrayTypeFunctionMapping, _boundary_to_dask, _pad_array,
-                             _validate_boundary, cuda_args, ngjit, not_implemented_func)
+                             _validate_boundary, _validate_raster, _validate_scalar,
+                             cuda_args, ngjit, not_implemented_func)
 from xrspatial.dataset_support import supports_dataset
 
 # TODO: Make convolution more generic with numba first-class functions.
@@ -278,6 +279,8 @@ def mean(agg, passes=1, excludes=[np.nan], name='mean', boundary='nan'):
         Dimensions without coordinates: dim_0, dim_1
     """
 
+    _validate_raster(agg, func_name='mean', name='agg')
+    _validate_scalar(passes, func_name='mean', name='passes', dtype=int, min_val=1)
     _validate_boundary(boundary)
     out = agg.data.astype(float)
     for i in range(passes):
@@ -507,12 +510,7 @@ def apply(raster, kernel, func=_calc_mean, name='focal_apply', boundary='nan'):
            [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_raster(raster, func_name='apply', name='raster')
 
     # Validate the kernel
     kernel = custom_kernel(kernel)
@@ -957,12 +955,7 @@ def focal_stats(agg,
           * 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_raster(agg, func_name='focal_stats', name='agg')
 
     # Validate the kernel
     kernel = custom_kernel(kernel)
@@ -1207,12 +1200,7 @@ def hotspots(raster, kernel, boundary='nan'):
         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")
+    _validate_raster(raster, func_name='hotspots', name='raster')
 
     _validate_boundary(boundary)
 

xrspatial/hillshade.py

@@ -14,6 +14,7 @@
 
 from .gpu_rtx import has_rtx
 from .utils import (_boundary_to_dask, _pad_array, _validate_boundary,
+                    _validate_raster, _validate_scalar,
                     calc_cuda_dims, get_dataarray_resolution,
                     has_cuda_and_cupy, is_cupy_array, is_cupy_backed)
 from .dataset_support import supports_dataset
@@ -217,6 +218,10 @@ def hillshade(agg: xr.DataArray,
         >>> hillshade_agg = hillshade(raster)
     """
 
+    _validate_raster(agg, func_name='hillshade', name='agg')
+    _validate_scalar(azimuth, func_name='hillshade', name='azimuth', min_val=0, max_val=360)
+    _validate_scalar(angle_altitude, func_name='hillshade', name='angle_altitude', min_val=0, max_val=90)
+
     if shadows and not has_rtx():
         raise RuntimeError(
             "Can only calculate shadows if cupy and rtxpy are available")

xrspatial/multispectral.py

@@ -9,8 +9,8 @@
 from numba import cuda
 from xarray import DataArray
 
-from xrspatial.utils import (ArrayTypeFunctionMapping, cuda_args, ngjit, not_implemented_func,
-                             validate_arrays)
+from xrspatial.utils import (ArrayTypeFunctionMapping, _validate_raster, cuda_args, ngjit,
+                             not_implemented_func, validate_arrays)
 from xrspatial.dataset_support import supports_dataset_bands
 
 # 3rd-party
@@ -153,6 +153,10 @@ def arvi(nir_agg: xr.DataArray,
          [ 0.02488688  0.00816024  0.00068681  0.02650602]]
     """
 
+    _validate_raster(nir_agg, func_name='arvi', name='nir_agg')
+    _validate_raster(red_agg, func_name='arvi', name='red_agg')
+    _validate_raster(blue_agg, func_name='arvi', name='blue_agg')
+
     validate_arrays(red_agg, nir_agg, blue_agg)
 
     mapper = ArrayTypeFunctionMapping(numpy_func=_arvi_cpu,
@@ -312,6 +316,10 @@ def evi(nir_agg: xr.DataArray,
          [-8.53211    5.486726   0.8394608  3.5043988]]
     """
 
+    _validate_raster(nir_agg, func_name='evi', name='nir_agg')
+    _validate_raster(red_agg, func_name='evi', name='red_agg')
+    _validate_raster(blue_agg, func_name='evi', name='blue_agg')
+
     if not red_agg.shape == nir_agg.shape == blue_agg.shape:
         raise ValueError("input layers expected to have equal shapes")
 
@@ -456,6 +464,9 @@ def gci(nir_agg: xr.DataArray,
          [0.34822243 0.28270411 0.29641694 0.359375  ]]
     """
 
+    _validate_raster(nir_agg, func_name='gci', name='nir_agg')
+    _validate_raster(green_agg, func_name='gci', name='green_agg')
+
     validate_arrays(nir_agg, green_agg)
 
     mapper = ArrayTypeFunctionMapping(numpy_func=_gci_cpu,
@@ -540,6 +551,9 @@ def nbr(nir_agg: xr.DataArray,
          [-0.10823033 -0.14486392 -0.12981689 -0.12121212]]
     """
 
+    _validate_raster(nir_agg, func_name='nbr', name='nir_agg')
+    _validate_raster(swir2_agg, func_name='nbr', name='swir2_agg')
+
     validate_arrays(nir_agg, swir2_agg)
 
     mapper = ArrayTypeFunctionMapping(
@@ -631,6 +645,9 @@ def nbr2(swir1_agg: xr.DataArray,
          [0.07218576 0.06857143 0.067659   0.07520281]]
     """
 
+    _validate_raster(swir1_agg, func_name='nbr2', name='swir1_agg')
+    _validate_raster(swir2_agg, func_name='nbr2', name='swir2_agg')
+
     validate_arrays(swir1_agg, swir2_agg)
 
     mapper = ArrayTypeFunctionMapping(
@@ -715,6 +732,9 @@ def ndvi(nir_agg: xr.DataArray,
          [0.06709956 0.04431737 0.04496226 0.07792632]]
     """
 
+    _validate_raster(nir_agg, func_name='ndvi', name='nir_agg')
+    _validate_raster(red_agg, func_name='ndvi', name='red_agg')
+
     validate_arrays(nir_agg, red_agg)
 
     mapper = ArrayTypeFunctionMapping(
@@ -804,6 +824,9 @@ def ndmi(nir_agg: xr.DataArray,
          [-0.17901748 -0.21133603 -0.19575651 -0.19464068]]
     """
 
+    _validate_raster(nir_agg, func_name='ndmi', name='nir_agg')
+    _validate_raster(swir1_agg, func_name='ndmi', name='swir1_agg')
+
     validate_arrays(nir_agg, swir1_agg)
 
     mapper = ArrayTypeFunctionMapping(
@@ -994,6 +1017,9 @@ def savi(nir_agg: xr.DataArray,
          [0.03353769 0.02215077 0.02247375 0.03895046]]
     """
 
+    _validate_raster(nir_agg, func_name='savi', name='nir_agg')
+    _validate_raster(red_agg, func_name='savi', name='red_agg')
+
     validate_arrays(red_agg, nir_agg)
 
     if not -1.0 <= soil_factor <= 1.0:
@@ -1138,6 +1164,10 @@ def sipi(nir_agg: xr.DataArray,
          [1.2903225 1.6451613 1.9708029 1.3556485]]
     """
 
+    _validate_raster(nir_agg, func_name='sipi', name='nir_agg')
+    _validate_raster(red_agg, func_name='sipi', name='red_agg')
+    _validate_raster(blue_agg, func_name='sipi', name='blue_agg')
+
     validate_arrays(red_agg, nir_agg, blue_agg)
 
     mapper = ArrayTypeFunctionMapping(numpy_func=_sipi_cpu,
@@ -1314,6 +1344,10 @@ def ebbi(red_agg: xr.DataArray,
             * lon      (lon) float64 0.0 1.0 2.0 3.0
     """
 
+    _validate_raster(red_agg, func_name='ebbi', name='red_agg')
+    _validate_raster(swir_agg, func_name='ebbi', name='swir_agg')
+    _validate_raster(tir_agg, func_name='ebbi', name='tir_agg')
+
     validate_arrays(red_agg, swir_agg, tir_agg)
 
     mapper = ArrayTypeFunctionMapping(numpy_func=_ebbi_cpu,
@@ -1515,6 +1549,10 @@ def true_color(r, g, b, nodata=1, c=10.0, th=0.125, name='true_color'):
         >>> true_color_img.plot.imshow()
     """
 
+    _validate_raster(r, func_name='true_color', name='r')
+    _validate_raster(g, func_name='true_color', name='g')
+    _validate_raster(b, func_name='true_color', name='b')
+
     mapper = ArrayTypeFunctionMapping(
         numpy_func=_true_color_numpy,
         dask_func=_true_color_dask,