Adding max distance to plotting script (#33)

CHANGELOG.rst

@@ -8,6 +8,7 @@ unreleased
 * Add zenodo release information (:pr:`24`) `Hauke Schulz`_.
 * Fix warning about seeting copy of a slice of a DataFrame (:pr:`28`) `Marius Rixen`_.
 * Update infrastructure to use pyproject.toml with pdm and ruff (:pr:`29`) `Hauke Schulz`_.
+* Added a function computing the distance tarvelled on the ascending branch and adding it to the trajectory plot. (:pr:`33`) `Marius Winkler`_.
 
 
 0.0.5 (2023-10-19)

pyproject.toml

@@ -22,6 +22,7 @@ dependencies = [
     "regex",
     "tqdm",
     "xarray",
+    "geopy",
 ]
 requires-python = ">=3.9"
 readme = "README.md"

scripts/plot_radiosonde_data.py

@@ -8,7 +8,7 @@
 @author: laura
 
 run the script:
-    python make_plots_radiosonde.py -i inputfilename -o outputdirectory
+    python plot_radiosonde_data.py -i inputfilename -o outputdirectory
 for level 1 radiosonde files. The inputfilename should include the direction (ascent or descent). If not, no SkewT diagram is made.
 
 This script produces 3 plots:
@@ -27,6 +27,7 @@
 import numpy as np
 import pint_xarray  # noqa: F401
 import xarray as xr
+from geopy.distance import geodesic
 from metpy.plots import Hodograph, SkewT
 from metpy.units import units
 
@@ -113,6 +114,36 @@ def data4skewt(direction, data):
     return p, T, Td, wind_speed, wind_dir, u, v
 
 
+def calculate_max_distance(ds, sounding):
+    latitudes = ds["lat"].values
+    longitudes = ds["lon"].values
+
+    # Ensure there are no NaN values by filtering them out
+    valid_indices = ~np.isnan(latitudes) & ~np.isnan(longitudes)
+    latitudes = latitudes[valid_indices]
+    longitudes = longitudes[valid_indices]
+
+    # Skip soundings with insufficient valid data
+    if len(latitudes) < 2 or len(longitudes) < 2:
+        print(f"Not enough valid data for sounding {sounding}. Skipping.")
+        return None
+
+    # Get the surface (initial) latitude and longitude
+    surface_lat = latitudes[0]
+    surface_lon = longitudes[0]
+
+    # Calculate the distance from the surface point to each other point
+    distances = np.array(
+        [
+            geodesic((surface_lat, surface_lon), (lat, lon)).kilometers
+            for lat, lon in zip(latitudes, longitudes)
+        ]
+    )
+
+    # Return the maximum distance traveled
+    return distances.max()
+
+
 def main():
     args = get_args()
 
@@ -194,13 +225,23 @@ def main():
     alt_max = Decimal(str(alt_max)).quantize(Decimal("0.1"), rounding=ROUND_HALF_UP)
     ax.text(
         0.02,
-        0.07,
+        0.14,
         f"max. height: {alt_max} km",
         transform=ax.transAxes,
         fontsize=11,
         verticalalignment="top",
         bbox=props,
     )
+    max_distance = calculate_max_distance(data, snd).round(1)
+    ax.text(
+        0.02,
+        0.07,
+        f"max. distance: {max_distance} km",
+        transform=ax.transAxes,
+        fontsize=11,
+        verticalalignment="top",
+        bbox=props,
+    )
     fig_track.savefig(
         f"{outdir}/Trajectories/track_{timestamp}_{direction}.png", bbox_inches="tight"
     )