diff --git a/xrspatial/classify.py b/xrspatial/classify.py
index d3261a5b9..b085a8488 100644
--- a/xrspatial/classify.py
+++ b/xrspatial/classify.py
@@ -842,12 +842,17 @@ def natural_breaks(agg: xr.DataArray,
         of values to be reclassified.
     k : int, default=5
         Number of classes to be produced.
-    num_sample : int, default=20000
+    num_sample : int or None, default=20000
         Number of sample data points used to fit the model.
         Natural Breaks (Jenks) classification is indeed O(n²) complexity,
         where n is the total number of data points, i.e: `agg.size`
         When n is large, we should fit the model on a small sub-sample
         of the data instead of using the whole dataset.
+        ``None`` means fit on all data instead of a sub-sample. That is
+        the full O(n²) case described above, so it may be slow and raises
+        ``MemoryError`` if the Jenks matrices would exceed half of the
+        available RAM. For dask the full sample is drawn lazily via
+        indexed access.
     name : str, default='natural_breaks'
         Name of output aggregate.