diff --git a/pyproject.toml b/pyproject.toml
index 7a011a0..3430df4 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -85,7 +85,7 @@ select = [
   "PLE",  # Pylint error https://github.com/charliermarsh/ruff#error-ple
 ]
 ignore = [
-"D100", "D101", "D104", "D105", "D106", "D107", "D203", "D213"
+"D100", "D101", "D104", "D105", "D106", "D107", "D203", "D213", "D413"
 ] # docstring style
 
 line-length = 88
diff --git a/sarxarray/stack.py b/sarxarray/stack.py
index b5f9299..09b671d 100644
--- a/sarxarray/stack.py
+++ b/sarxarray/stack.py
@@ -58,6 +58,9 @@ def point_selection(self, threshold, method="amplitude_dispersion", chunks=1000)
         """
         match method:
             case "amplitude_dispersion":
+                # Amplitude dispersion thresholding
+                # Note there can be NaN values in the amplitude dispersion
+                # However NaN values will not pass this threshold
                 mask = self._amp_disp() < threshold
             case _:
                 raise NotImplementedError
@@ -111,9 +114,14 @@ def _amp_disp(self, chunk_azimuth=500, chunk_range=500):
             {"azimuth": chunk_azimuth, "range": chunk_range, "time": -1}
         )
 
-        amplitude_dispersion = amplitude.std(axis=t_order) / (
-            amplitude.mean(axis=t_order) + np.finfo(amplitude.dtype).eps
-        )  # adding epsilon to avoid zero division
+        # Compoute amplitude dispersion
+        # By defalut, the mean and std function from Xarray will skip NaN values
+        # However, if there is NaN value in time series, we want to discard the pixel
+        # Therefore, we set skipna=False
+        # Adding epsilon to avoid zero division
+        amplitude_dispersion = amplitude.std(axis=t_order, skipna=False) / (
+            amplitude.mean(axis=t_order, skipna=False) + np.finfo(amplitude.dtype).eps
+        )
 
         return amplitude_dispersion
 
diff --git a/tests/test_stack.py b/tests/test_stack.py
index 0af56d7..272ab71 100644
--- a/tests/test_stack.py
+++ b/tests/test_stack.py
@@ -72,7 +72,7 @@ def test_stack_pointselection_some(self, synthetic_dataset):
         synthetic_dataset = synthetic_dataset.slcstack._get_phase()
         stm = synthetic_dataset.slcstack.point_selection(
             threshold=0.2, method="amplitude_dispersion"
-        )  # select all
+        )  # select some
         assert stm.space.shape[0] == 4
         assert set([k for k in stm.coords.keys()]).issubset(
             ["time", "azimuth", "range"]
@@ -81,6 +81,35 @@ def test_stack_pointselection_some(self, synthetic_dataset):
             ["complex", "amplitude", "phase"]
         )
 
+    def test_stack_pointselection_nan(self):
+        data = np.random.rand(10, 10, 10) + 1j * np.random.rand(10, 10, 10)
+        data[0, 0, 0] = np.nan  # introduce a nan
+        ds_nan = xr.Dataset(
+            {
+                "complex": (
+                    ("azimuth", "range", "time"),
+                    data,
+                )
+            },
+            coords={
+                "azimuth": np.arange(600, 610, 1, dtype=int),
+                "range": np.arange(1400, 1410, 1, dtype=int),
+                "time": np.arange(1, 11, 1, dtype=int),
+            },
+        )
+        ds_nan = ds_nan.slcstack._get_amplitude()
+        ds_nan = ds_nan.slcstack._get_phase()
+        stm = ds_nan.slcstack.point_selection(
+            threshold=100, method="amplitude_dispersion"
+        )
+        assert stm.space.shape[0] == 99  # only the nan is removed
+        assert set(list(stm.coords.keys())).issubset(
+            ["time", "azimuth", "range"]
+        )
+        assert set(list(stm.data_vars.keys())).issubset(
+            ["complex", "amplitude", "phase"]
+        )
+
 
 class TestStackMultiLook:
     def test_stack_multi_look_mean(self, synthetic_dataset):