Repository: makepath/xarray-spatial

Issue:
unexpected behaviour for `natural_breaks` with sampling
Thanks a lot for this fantastic library!   

**Describe the bug**
I have observed the following behaviour, which I think is unexpected, when computing Jenks natural breaks on a sub-sample of points. 
All points that are larger than the largest point in the sub-sample seems to be classified as NaN. See code snippet below: if the sub-sample includes the largest (finite) value in the dataset, then all points are properly classified. However, if the sub-sample only covers a part of the data range, one is left with some of the valid points being misclassified.

**Expected behavior**
I would expect that even if the breaks are computed on a sub-sample, all finite points (including points that fall outside the range of the sub-sample) should be classified. 

**Screenshots**
```python
import numpy as np
import xarray as xr

from xrspatial.classify import natural_breaks

raster = xr.DataArray(
    np.array([
        [1.,  2.,  3.,  4., np.nan],
        [5.,  6.,  7.,  8.,  9.],
        [10., 11., 12., 13., 14.],
        [15., 16., 17., 18., np.inf],
    ])
)

#  Here all finite values are properly labeled
agg = natural_breaks(raster, k=5, num_sample=25)
print(agg.values)
# [[ 0.  0.  0.  1. nan]
#  [ 1.  1.  2.  2.  2.]
#  [ 2.  3.  3.  3.  3.]
#  [ 4.  4.  4.  4. nan]]

# Here the value in (3,3) is classified as "NaN" - I would have expected "4"
agg = natural_breaks(raster, k=5, num_sample=8)
print(agg.values)
# [[ 0.  0.  0.  0. nan]
#  [ 0.  1.  1.  1.  2.]
#  [ 2.  3.  3.  3.  3.]
#  [ 4.  4.  4. nan nan]]
```

**Desktop (please complete the following information):**
 - OS X
 - Version 0.3.2

**Additional context**
At first I thought that this issue could be fixed by changing the upper-bounds for the classification:
https://github.com/makepath/xarray-spatial/blob/1f2282afdd3eeac201b61be491d1a09a4a0ce2b8/xrspatial/classify.py#L530
to:
```python
kclass[-1] = np.inf
```
However, this would lead to also `np.inf` in the input data to be assigned to the class with label `4` - but maybe `np.inf`/`-np.inf` should also actually be classified?  


Discussion:

Patch:
diff --git a/benchmarks/benchmarks/aspect.py b/benchmarks/benchmarks/aspect.py
index 7147671..3698090 100644
--- a/benchmarks/benchmarks/aspect.py
+++ b/benchmarks/benchmarks/aspect.py
@@ -1,14 +1,10 @@
 from xrspatial import aspect
-from .common import get_xr_dataarray
+from .common import Benchmarking
 
 
-class Aspect:
-    params = ([100, 300, 1000, 3000, 10000], ["numpy", "cupy"])
-    param_names = ("nx", "type")
-
-    def setup(self, nx, type):
-        ny = nx // 2
-        self.xr = get_xr_dataarray((ny, nx), type)
+class Aspect(Benchmarking):
+    def __init__(self):
+        super().__init__(func=aspect)
 
     def time_aspect(self, nx, type):
-        aspect(self.xr)
+        return self.time(nx, type)
diff --git a/benchmarks/benchmarks/classify.py b/benchmarks/benchmarks/classify.py
new file mode 100644
index 0000000..eb10cf1
--- /dev/null
+++ b/benchmarks/benchmarks/classify.py
@@ -0,0 +1,40 @@
+import numpy as np
+
+from xrspatial.classify import reclassify, equal_interval, quantile, natural_breaks
+from .common import get_xr_dataarray
+
+
+class Classify:
+    params = ([100, 300, 1000, 3000, 10000], [1, 10, 100], ["numpy", "cupy"])
+    param_names = ("nx", "nbins", "type")
+
+    def setup(self, nx, nbins, type):
+        ny = nx // 2
+        self.agg = get_xr_dataarray((ny, nx), type)
+        min_val = np.nanmin(self.agg.data)
+        max_val = np.nanmax(self.agg.data)
+        self.nbins = nbins
+        self.bins = np.linspace(min_val, max_val, self.nbins)
+        self.new_values = np.arange(nbins)
+
+
+class Reclassify(Classify):
+    def time_reclassify(self, nx, nbins, type):
+        reclassify(self.agg, self.bins, self.new_values)
+
+
+class Quantile(Classify):
+    def time_quantile(self, nx, nbins, type):
+        quantile(self.agg, k=self.nbins)
+
+
+class NaturalBreaks(Classify):
+    params = ([100, 300, 1000, 3000, 10000], [1, 10], ["numpy", "cupy"])
+
+    def time_natural_breaks(self, nx, nbins, type):
+        natural_breaks(self.agg, k=self.nbins)
+
+
+class EqualInterval(Classify):
+    def time_equal_interval(self, nx, nbins, type):
+        equal_interval(self.agg, k=self.nbins)
diff --git a/benchmarks/benchmarks/common.py b/benchmarks/benchmarks/common.py
index ac0b738..a125320 100644
--- a/benchmarks/benchmarks/common.py
+++ b/benchmarks/benchmarks/common.py
@@ -4,7 +4,9 @@ from xrspatial.gpu_rtx import has_rtx
 from xrspatial.utils import has_cuda, has_cupy
 
 
-def get_xr_dataarray(shape, type, different_each_call=False):
+def get_xr_dataarray(
+    shape, type, different_each_call=False, seed=71942, is_int=False, include_nan=False
+):
     # Gaussian bump with noise.
     #
     # Valid types are "numpy", "cupy" and "rtxpy". Using "numpy" will return
@@ -20,16 +22,25 @@ def get_xr_dataarray(shape, type, different_each_call=False):
     # optimistically fast benchmark times.
     ny, nx = shape
 
-    x = np.linspace(-1000, 1000, nx)
-    y = np.linspace(-800, 800, ny)
+    x = np.linspace(-180, 180, nx)
+    y = np.linspace(-90, 90, ny)
     x2, y2 = np.meshgrid(x, y)
-    z = 100.0*np.exp(-x2**2 / 5e5 - y2**2 / 2e5)
+    rng = np.random.default_rng(seed)
 
-    rng = np.random.default_rng(71942)
-    z += rng.normal(0.0, 2.0, (ny, nx))
+    if is_int:
+        z = rng.integers(-nx, nx, size=shape).astype(np.float32)
+    else:
+        z = 100.0*np.exp(-x2**2 / 5e5 - y2**2 / 2e5)
+        z += rng.normal(0.0, 2.0, (ny, nx))
 
     if different_each_call:
-        z[-1, -1] = np.random.default_rng().normal(0.0, 2.0)
+        if is_int:
+            z[-1, -1] = np.random.default_rng().integers(-nx, nx)
+        else:
+            z[-1, -1] = np.random.default_rng().normal(0.0, 2.0)
+
+    if include_nan:
+        z[0, 0] = np.nan
 
     if type == "numpy":
         pass
@@ -46,4 +57,20 @@ def get_xr_dataarray(shape, type, different_each_call=False):
     else:
         raise RuntimeError(f"Unrecognised type {type}")
 
-    return xr.DataArray(z, coords=dict(x=x, y=y), dims=["y", "x"])
+    return xr.DataArray(z, coords=dict(y=y, x=x), dims=["y", "x"])
+
+
+class Benchmarking:
+    params = ([100, 300, 1000, 3000, 10000], ["numpy", "cupy"])
+    param_names = ("nx", "type")
+
+    def __init__(self, func=None):
+        self.func = func
+
+    def setup(self, nx, type):
+        ny = nx // 2
+        self.xr = get_xr_dataarray((ny, nx), type)
+
+    def time(self, nx, type):
+        if self.func is not None:
+            self.func(self.xr)
diff --git a/benchmarks/benchmarks/curvature.py b/benchmarks/benchmarks/curvature.py
new file mode 100644
index 0000000..003b359
--- /dev/null
+++ b/benchmarks/benchmarks/curvature.py
@@ -0,0 +1,10 @@
+from xrspatial import curvature
+from .common import Benchmarking
+
+
+class Curvature(Benchmarking):
+    def __init__(self):
+        super().__init__(func=curvature)
+
+    def time_curvature(self, nx, type):
+        return self.time(nx, type)
diff --git a/benchmarks/benchmarks/focal.py b/benchmarks/benchmarks/focal.py
new file mode 100644
index 0000000..3cdb224
--- /dev/null
+++ b/benchmarks/benchmarks/focal.py
@@ -0,0 +1,45 @@
+import numpy as np
+
+from xrspatial.convolution import custom_kernel
+from xrspatial.focal import mean, apply, hotspots, focal_stats
+from .common import get_xr_dataarray
+
+
+class Focal:
+    params = ([100, 300, 1000, 3000], [(5, 5), (25, 25)], ["numpy", "cupy"])
+    param_names = ("nx", "kernelsize", "type")
+
+    def setup(self, nx, kernelsize, type):
+        ny = nx // 2
+        self.agg = get_xr_dataarray((ny, nx), type)
+        kernel_w, kernel_h = kernelsize
+        self.kernel = custom_kernel(np.ones((kernel_h, kernel_w)))
+
+
+class FocalApply(Focal):
+    def time_apply(self, nx, kernelsize, type):
+        apply(self.agg, self.kernel)
+
+
+class FocalHotspots(Focal):
+    def time_hotspots(self, nx, kernelsize, type):
+        hotspots(self.agg, self.kernel)
+
+
+class FocalStats(Focal):
+    params = ([100, 300, 1000, 3000], [(5, 5), (15, 15)], ["numpy", "cupy"])
+
+    def time_focal_stats(self, nx, kernelsize, type):
+        focal_stats(self.agg, self.kernel)
+
+
+class FocalMean:
+    params = ([100, 300, 1000, 3000, 10000], [1, 10], ["numpy", "cupy"])
+    param_names = ("nx", "passes", "type")
+
+    def setup(self, nx, passes, type):
+        ny = nx // 2
+        self.agg = get_xr_dataarray((ny, nx), type)
+
+    def time_mean(self, nx, passes, type):
+        mean(self.agg, passes)
diff --git a/benchmarks/benchmarks/multispectral.py b/benchmarks/benchmarks/multispectral.py
new file mode 100644
index 0000000..fd86e01
--- /dev/null
+++ b/benchmarks/benchmarks/multispectral.py
@@ -0,0 +1,67 @@
+from xrspatial.multispectral import arvi, evi, gci, nbr, nbr2, ndvi, ndmi, savi, sipi, ebbi
+from .common import get_xr_dataarray, Benchmarking
+
+
+class Multispectral(Benchmarking):
+    def __init__(self):
+        super().__init__()
+
+    def setup(self, nx, type):
+        ny = nx // 2
+        self.red = get_xr_dataarray((ny, nx), type, seed=100)
+        self.green = get_xr_dataarray((ny, nx), type, seed=200)
+        self.blue = get_xr_dataarray((ny, nx), type, seed=300)
+        self.nir = get_xr_dataarray((ny, nx), type, seed=400)
+        self.swir1 = get_xr_dataarray((ny, nx), type, seed=500)
+        self.swir2 = get_xr_dataarray((ny, nx), type, seed=600)
+        self.tir = get_xr_dataarray((ny, nx), type, seed=700)
+
+
+class Arvi(Multispectral):
+    def time_arvi(self, nx, type):
+        arvi(self.nir, self.red, self.blue)
+
+
+class Evi(Multispectral):
+    def time_evi(self, nx, type):
+        evi(self.nir, self.red, self.blue)
+
+
+class Gci(Multispectral):
+    def time_gci(self, nx, type):
+        gci(self.nir, self.green)
+
+
+class Nbr(Multispectral):
+    def time_nbr(self, nx, type):
+        nbr(self.nir, self.swir2)
+
+
+class Nbr2(Multispectral):
+    def time_nbr2(self, nx, type):
+        nbr2(self.swir1, self.swir2)
+
+
+class Ndvi(Multispectral):
+    def time_ndvi(self, nx, type):
+        ndvi(self.nir, self.red)
+
+
+class Ndmi(Multispectral):
+    def time_ndmi(self, nx, type):
+        ndmi(self.nir, self.swir1)
+
+
+class Savi(Multispectral):
+    def time_savi(self, nx, type):
+        savi(self.nir, self.red)
+
+
+class Sipi(Multispectral):
+    def time_sipi(self, nx, type):
+        sipi(self.nir, self.red, self.blue)
+
+
+class Ebbi(Multispectral):
+    def time_ebbi(self, nx, type):
+        ebbi(self.red, self.swir1, self.tir)
diff --git a/benchmarks/benchmarks/pathfinding.py b/benchmarks/benchmarks/pathfinding.py
new file mode 100644
index 0000000..1dc84bc
--- /dev/null
+++ b/benchmarks/benchmarks/pathfinding.py
@@ -0,0 +1,20 @@
+from xrspatial.pathfinding import a_star_search
+from .common import get_xr_dataarray
+
+
+class AStarSearch:
+    params = ([10, 100, 300], [4, 8], ["numpy"])
+    param_names = ("nx", "connectivity", "type")
+
+    def setup(self, nx, connectivity, type):
+        ny = nx // 2
+        self.agg = get_xr_dataarray((ny, nx), type)
+        self.start = self.agg.y[0], self.agg.x[0]
+        self.goal = self.agg.y[-1], self.agg.x[-1]
+
+    def time_a_star_search(self, nx, connectivity, type):
+        a_star_search(
+            self.agg, self.start, self.goal,
+            connectivity=connectivity,
+            snap_start=True, snap_goal=True
+        )
diff --git a/benchmarks/benchmarks/perlin.py b/benchmarks/benchmarks/perlin.py
new file mode 100644
index 0000000..ac2b063
--- /dev/null
+++ b/benchmarks/benchmarks/perlin.py
@@ -0,0 +1,11 @@
+from xrspatial.perlin import perlin
+
+from .common import Benchmarking
+
+
+class Perlin(Benchmarking):
+    def __init__(self):
+        super().__init__(func=perlin)
+
+    def time_perlin(self, nx, type):
+        return self.time(nx, type)
diff --git a/benchmarks/benchmarks/proximity.py b/benchmarks/benchmarks/proximity.py
new file mode 100644
index 0000000..f449845
--- /dev/null
+++ b/benchmarks/benchmarks/proximity.py
@@ -0,0 +1,35 @@
+import numpy as np
+
+from xrspatial.proximity import proximity, allocation, direction
+from .common import get_xr_dataarray
+
+
+class Base:
+    params = (
+        [100, 1000],
+        [1, 10, 100],
+        ["EUCLIDEAN", "GREAT_CIRCLE", "MANHATTAN"],
+        ["numpy"]
+    )
+    param_names = ("nx", "n_target_values", "distance_metric", "type")
+
+    def setup(self, nx, n_target_values, distance_metric, type):
+        ny = nx // 2
+        self.agg = get_xr_dataarray((ny, nx), type, is_int=True)
+        unique_values = np.unique(self.agg.data)
+        self.target_values = unique_values[:n_target_values]
+
+
+class Proximity(Base):
+    def time_proximity(self, nx, n_target_values, distance_metric, type):
+        proximity(self.agg, target_values=self.target_values, distance_metric=distance_metric)
+
+
+class Allocation(Base):
+    def time_allocation(self, nx, n_target_values, distance_metric, type):
+        allocation(self.agg, target_values=self.target_values, distance_metric=distance_metric)
+
+
+class Direction(Base):
+    def time_direction(self, nx, n_target_values, distance_metric, type):
+        direction(self.agg, target_values=self.target_values, distance_metric=distance_metric)
diff --git a/benchmarks/benchmarks/slope.py b/benchmarks/benchmarks/slope.py
index 6a1f528..101c3e8 100644
--- a/benchmarks/benchmarks/slope.py
+++ b/benchmarks/benchmarks/slope.py
@@ -1,14 +1,10 @@
 from xrspatial import slope
-from .common import get_xr_dataarray
+from .common import Benchmarking
 
 
-class Slope:
-    params = ([100, 300, 1000, 3000, 10000], ["numpy", "cupy"])
-    param_names = ("nx", "type")
-
-    def setup(self, nx, type):
-        ny = nx // 2
-        self.xr = get_xr_dataarray((ny, nx), type)
+class Slope(Benchmarking):
+    def __init__(self):
+        super().__init__(func=slope)
 
     def time_slope(self, nx, type):
-        slope(self.xr)
+        return self.time(nx, type)
diff --git a/benchmarks/benchmarks/terrain.py b/benchmarks/benchmarks/terrain.py
new file mode 100644
index 0000000..9ae956d
--- /dev/null
+++ b/benchmarks/benchmarks/terrain.py
@@ -0,0 +1,11 @@
+from xrspatial.terrain import generate_terrain
+
+from .common import Benchmarking
+
+
+class GenerateTerrain(Benchmarking):
+    def __init__(self):
+        super().__init__(func=generate_terrain)
+
+    def time_generate_terrain(self, nx, type):
+        return self.time(nx, type)
diff --git a/xrspatial/classify.py b/xrspatial/classify.py
index 7a09f44..a2492f6 100755
--- a/xrspatial/classify.py
+++ b/xrspatial/classify.py
@@ -16,8 +16,6 @@ from datashader.colors import rgb
 import numba as nb
 import dask.array as da
 
-from numpy.random import RandomState
-
 from xrspatial.utils import cuda_args
 from xrspatial.utils import ngjit
 from xrspatial.utils import ArrayTypeFunctionMapping
@@ -78,11 +76,13 @@ def _cpu_bin(data, bins, new_values):
                 # first bin
                 if b == 0:
                     if val <= bins[b]:
-                        val_bin = b
+                        if np.isfinite(val):
+                            val_bin = b
                         break
                 else:
                     if val > bins[b - 1] and val <= bins[b]:
-                        val_bin = b
+                        if np.isfinite(val):
+                            val_bin = b
                         break
 
             if val_bin > -1:
@@ -144,6 +144,10 @@ def _run_gpu_bin(data, bins, new_values, out):
 
 
 def _run_cupy_bin(data, bins, new_values):
+    # replace inf by nan to avoid classify these values as we want to treat them as outliers
+    data = cupy.where(data == cupy.inf, cupy.nan, data)
+    data = cupy.where(data == -cupy.inf, cupy.nan, data)
+
     bins_cupy = cupy.asarray(bins)
     new_values_cupy = cupy.asarray(new_values)
     out = cupy.empty(data.shape, dtype='f4')
@@ -542,12 +546,17 @@ def _run_jenks(data, n_classes, module):
 def _run_natural_break(agg, num_sample, k, module):
     data = agg.data
     num_data = data.size
+    max_data = module.max(data[module.isfinite(data)])
+    if module == cupy:
+        max_data = max_data.get()
 
     if num_sample is not None and num_sample < num_data:
         # randomly select sample from the whole dataset
         # create a pseudo random number generator
-        generator = RandomState(1234567890)
-        idx = module.linspace(
+        # Note: cupy and nupy generate different random numbers
+        # use numpy.random to ensure the same result
+        generator = np.random.RandomState(1234567890)
+        idx = np.linspace(
             0, data.size, data.size, endpoint=False, dtype=module.uint32
         )
         generator.shuffle(idx)
@@ -586,6 +595,7 @@ def _run_natural_break(agg, num_sample, k, module):
     else:
         centroids = _run_jenks(sample_data, k, module)
         bins = module.array(centroids[1:])
+        bins[-1] = max_data
 
     out = _bin(agg, bins, module.arange(uvk))
     return out
@@ -701,6 +711,8 @@ def _run_equal_interval(agg, k, module):
         inf = np.inf
 
     data = module.where(data == inf, nan, data)
+    data = module.where(data == -inf, nan, data)
+
     max_data = module.nanmax(data)
     min_data = module.nanmin(data)
     if module == cupy:


Test patch:
diff --git a/xrspatial/tests/test_classify.py b/xrspatial/tests/test_classify.py
index 55fc2bb..073ca66 100755
--- a/xrspatial/tests/test_classify.py
+++ b/xrspatial/tests/test_classify.py
@@ -15,7 +15,7 @@ from xrspatial.tests.general_checks import create_test_raster
 
 def input_data(backend='numpy'):
     elevation = np.array([
-        [1.,  2.,  3.,  4., np.nan],
+        [-np.inf,  2.,  3.,  4., np.nan],
         [5.,  6.,  7.,  8.,  9.],
         [10., 11., 12., 13., 14.],
         [15., 16., 17., 18., np.inf],
@@ -29,10 +29,10 @@ def result_reclassify():
     bins = [10, 15, np.inf]
     new_values = [1, 2, 3]
     expected_result = np.asarray([
-        [1., 1., 1., 1., np.nan],
+        [np.nan, 1., 1., 1., np.nan],
         [1., 1., 1., 1., 1.],
         [1., 2., 2., 2., 2.],
-        [2., 3., 3., 3., 3.]
+        [2., 3., 3., 3., np.nan]
     ], dtype=np.float32)
     return bins, new_values, expected_result
 
@@ -71,8 +71,8 @@ def test_reclassify_dask_cupy(result_reclassify):
 def result_quantile():
     k = 5
     expected_result = np.asarray([
-        [0., 0., 0., 0., np.nan],
-        [1., 1., 1., 2., 2.],
+        [np.nan, 0., 0., 0., np.nan],
+        [0., 1., 1., 1., 2.],
         [2., 2., 3., 3., 3.],
         [4., 4., 4., 4., np.nan]
     ], dtype=np.float32)
@@ -117,14 +117,27 @@ def test_quantile_cupy(result_quantile):
 def result_natural_breaks():
     k = 5
     expected_result = np.asarray([
-        [0., 0., 0., 1., np.nan],
-        [1., 1., 2., 2., 2.],
+        [np.nan, 0., 0., 0., np.nan],
+        [1., 1., 1., 2., 2.],
         [2., 3., 3., 3., 3.],
         [4., 4., 4., 4., np.nan]
     ], dtype=np.float32)
     return k, expected_result
 
 
+@pytest.fixture
+def result_natural_breaks_num_sample():
+    k = 5
+    num_sample = 8
+    expected_result = np.asarray([
+        [np.nan, 0., 0., 0., np.nan],
+        [0., 1., 1., 1., 2.],
+        [2., 3., 3., 3., 3.],
+        [4., 4., 4., 4., np.nan]
+    ], dtype=np.float32)
+    return k, num_sample, expected_result
+
+
 def test_natural_breaks_numpy(result_natural_breaks):
     numpy_agg = input_data()
     k, expected_result = result_natural_breaks
@@ -132,6 +145,13 @@ def test_natural_breaks_numpy(result_natural_breaks):
     general_output_checks(numpy_agg, numpy_natural_breaks, expected_result, verify_dtype=True)
 
 
+def test_natural_breaks_numpy_num_sample(result_natural_breaks_num_sample):
+    numpy_agg = input_data()
+    k, num_sample, expected_result = result_natural_breaks_num_sample
+    numpy_natural_breaks = natural_breaks(numpy_agg, k=k, num_sample=num_sample)
+    general_output_checks(numpy_agg, numpy_natural_breaks, expected_result, verify_dtype=True)
+
+
 def test_natural_breaks_cpu_deterministic():
     results = []
     elevation = np.arange(100).reshape(10, 10)
@@ -164,12 +184,20 @@ def test_natural_breaks_cupy(result_natural_breaks):
     general_output_checks(cupy_agg, cupy_natural_breaks, expected_result, verify_dtype=True)
 
 
+@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available")
+def test_natural_breaks_cupy_num_sample(result_natural_breaks_num_sample):
+    cupy_agg = input_data('cupy')
+    k, num_sample, expected_result = result_natural_breaks_num_sample
+    cupy_natural_breaks = natural_breaks(cupy_agg, k=k, num_sample=num_sample)
+    general_output_checks(cupy_agg, cupy_natural_breaks, expected_result, verify_dtype=True)
+
+
 @pytest.fixture
 def result_equal_interval():
     k = 3
     expected_result = np.asarray([
-        [0., 0., 0., 0., np.nan],
-        [0., 0., 1., 1., 1.],
+        [np.nan, 0., 0., 0., np.nan],
+        [0., 0., 0., 1., 1.],
         [1., 1., 1., 2., 2.],
         [2., 2., 2., 2., np.nan]
     ], dtype=np.float32)