From 44690e1b209fcfa5c3acbd2b3001573918f0cae1 Mon Sep 17 00:00:00 2001
From: William Kearney <william.kearney@uni-potsdam.de>
Date: Thu, 28 May 2026 09:54:02 +0200
Subject: [PATCH] PPS: Add plotdz method

This is slightly different from the MATLAB PPS.plotdz method in that
it doesn't plot the stream profile. It is tricky to forward some of
the necessary keyword arguments to StreamObject.plotdz and some to
plt.scatter. Use

...
p = tt3.PPS.from_nal(s, kp)
p.s.plotdz(dem, ax=ax)
p.plotdz(dem, ax=ax)
...

to plot both the StreamObject and the PPS.

Signed-off-by: William Kearney <william.kearney@uni-potsdam.de>
---
 src/topotoolbox/pps.py | 95 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 95 insertions(+)

diff --git a/src/topotoolbox/pps.py b/src/topotoolbox/pps.py
index 0017a58..ecbeafb 100644
--- a/src/topotoolbox/pps.py
+++ b/src/topotoolbox/pps.py
@@ -3,6 +3,7 @@
 
 from dataclasses import dataclass
 
+import matplotlib.pyplot as plt
 import numpy as np
 import numpy.typing as npt
 
@@ -52,3 +53,97 @@ def from_nal(cls, stream: StreamObject, nal: npt.NDArray[np.bool]):
         """Construct a PPS from a logical node attribute list
         """
         return cls(stream, np.flatnonzero(nal))
+
+    def plotdz(self, z, distance=None,
+               ax=None,
+               dunit: str = 'm', doffset: float = 0,
+               scalex=True, scaley=True,
+               **kwargs):
+        """Plot upstream distances of points in a PPS against a covariate
+
+        Note that this will only plot the points of the PPS. To
+        additionally plot the stream profiles, use `p.s.plotdz(z,
+        distance)` separately.
+
+        Parameters
+        ----------
+        z: GridObject, np.ndarray
+          The node attribute list that will be plotted on the y axis
+
+        distance: GridObject, np.ndarray
+          The node attribute list that will be plotted on the x axis
+
+        ax: matplotlib.axes.Axes, optional
+            The axes in which to plot the StreamObject. If no axes are
+            given, the current axes are used.
+
+        dunit: str, optional
+            The unit to plot the upstream distance. Should be either
+            'm' for meters or 'km' for kilometers.
+
+        doffset: float, optional
+            An offset to be applied to the upstream distance.
+            `doffset` should be in the units specified by `dunit`.
+
+        scalex: bool, optional
+            Autoscale the x-axis limits. Defaults to `True`. If the
+            x-axis limits have been set manually with `set_xlim`, the
+            autoscaling will not be applied regardless of the value of
+            `scalex`.
+
+        scaley: bool, optional
+            Autoscale the y-axis limits. Defaults to `True`. If the
+            y-axis limits have been set manually with `set_ylim`, the
+            autoscaling will not be applied regardless of the value of
+            `scaley`.
+
+        **kwargs
+            Additional keyword arguments are forwarded to `scatter.`
+
+        Returns
+        -------
+        matplotlib.axes.Axes
+            The axes into which the plot as been added
+
+        Raises
+        ------
+        ValueError
+            If `dunit` is not one of 'm' or 'km'.
+
+        Example
+        -------
+        .. plot ::
+
+           >>> import topotoolbox
+           >>> import matplotlib.pyplot as plt
+           >>> dem = topotoolbox.load_dem('bigtujunga')
+           >>> fd = topotoolbox.FlowObject(dem)
+           >>> s = topotoolbox.StreamObject(fd, threshold=1000)
+           >>> s = s.klargestconncomps(1)
+           >>> kp = s.knickpointfinder(dem, tolerance=20.0)
+           >>> p = topotoolbox.PPS.from_nal(s, kp)
+           >>> fig, ax = plt.subplots()
+           >>> _ = s.plotdz(dem, ax=ax)
+           >>> _ = p.plotdz(dem, ax=ax, color='k', zorder=2)
+           >>> plt.show()
+        """
+        if ax is None:
+            ax = plt.gca()
+
+        if distance is None:
+            distance = self.s.upstream_distance()
+
+        dist = self.s.ezgetnal(distance)
+
+        if dunit == 'km':
+            dist /= 1000
+        elif dunit != 'm':
+            raise ValueError("dunit must be one of 'm' or 'km'")
+
+        dist += doffset
+
+        z = self.s.ezgetnal(z)
+        ax.scatter(dist[self.pp], z[self.pp], **kwargs)
+        ax.autoscale_view(scalex=scalex, scaley=scaley)
+
+        return ax