Merge pull request #39 from MeteoSwiss-APN/dev-Peiderk

Fix issue 38

Author: Peiderk

README.md

@@ -32,10 +32,9 @@ instructions here below to set up a conda environment and test multiple use-case
 
 ### Create environment
 
-Create an environment and install the package dependencies with the  script `setup_env.sh` provided in the `tools` directory.
-Check available options with
+Create an environment and install the package dependencies with conda
 
-    tools/setup_env.sh -h
+     conda env create -n pytrajplot -f requirements/environment.yml
 
 We distinguish pinned installations based on exported (reproducible) environments,
 saved in `requirements/environment.yml`,
@@ -72,7 +71,7 @@ installed by conda and should not be modified by pip, use the `--no-deps` flag.
 
 For development, install the package in editable mode:
 
-    pip install --editable --no-deps .
+    pip install --no-deps --editable .
 
 *Warning:* Make sure you use the right pip, i.e. the one from the installed conda environment (`which pip` should point to something like `path/to/miniconda/envs/<package_env_name>/bin/pip`).
 

src/pytrajplot/parse_data.py

@@ -360,34 +360,34 @@ def read_trajectory(trajectory_file_path, start_df, plot_info_dict):
     # clean up the df in case its generated from a COSMO trajectory file
     if case == "COSMO":
         # at missing data points (i.e. if trajectory leaves computational domain), the z & hsurf values default to -999
-        traj_df.loc[(traj_df["z"] < 0), "z"] = np.NaN
-        traj_df.loc[(traj_df["hsurf"] < 0), "hsurf"] = np.NaN
+        traj_df.loc[(traj_df["z"] < 0), "z"] = np.nan
+        traj_df.loc[(traj_df["hsurf"] < 0), "hsurf"] = np.nan
         traj_df.dropna(
             subset=["lon"], inplace=True
         )  # remove rows containing only the origin/z_type
         traj_df.reset_index(inplace=True)
 
     # clean up the df in case its generated from a HRES trajectory file
     if case == "HRES":
-        traj_df.loc[(traj_df["lon"] == -999.00), "lon"] = np.NaN
-        traj_df.loc[(traj_df["lat"] == -999.00), "lat"] = np.NaN
-        traj_df.loc[(traj_df["z"] == -999), "z"] = np.NaN
+        traj_df.loc[(traj_df["lon"] == -999.00), "lon"] = np.nan
+        traj_df.loc[(traj_df["lat"] == -999.00), "lat"] = np.nan
+        traj_df.loc[(traj_df["z"] == -999), "z"] = np.nan
 
     # add various (empty) keys to the trajectory dataframe
     traj_df["z_type"] = None
     traj_df["origin"] = None
     traj_df["side_traj"] = None
-    traj_df["start_altitude"] = np.NaN
-    traj_df["lon_precise"] = np.NaN
-    traj_df["lat_precise"] = np.NaN
+    traj_df["start_altitude"] = np.nan
+    traj_df["lon_precise"] = np.nan
+    traj_df["lat_precise"] = np.nan
     traj_df["altitude_levels"] = None
     traj_df["#trajectories"] = number_of_trajectories
     traj_df["block_length"] = number_of_times
     traj_df["trajectory_direction"] = trajectory_file_path[
         -1:
     ]  # last letter of key = F/B --> direction of trajectory
-    traj_df["subplot_index"] = np.NaN
-    traj_df["max_start_altitude"] = np.NaN
+    traj_df["subplot_index"] = np.nan
+    traj_df["max_start_altitude"] = np.nan
 
     # add information to newly created (empty) keys in trajectory dataframe
     # basically, at this point the information from the start_df, gets merged into the trajectory dataframe
@@ -443,8 +443,8 @@ def read_trajectory(trajectory_file_path, start_df, plot_info_dict):
         reference_time=reference_time,
     )
 
-    # change the lon/lat values where the trajectory leaves the domain from their computational domain-boundary values to np.NaN.
-    traj_df.loc[np.isnan(traj_df["z"]), ["lon", "lat"]] = np.NaN
+    # change the lon/lat values where the trajectory leaves the domain from their computational domain-boundary values to np.nan.
+    traj_df.loc[np.isnan(traj_df["z"]), ["lon", "lat"]] = np.nan
 
     return traj_df
 

src/pytrajplot/plotting/analyse_trajectories.py

@@ -9,7 +9,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
-from numpy.lib.function_base import mean
+from numpy import mean
 
 
 # FUNCTIONS
@@ -306,7 +306,8 @@ def _create_plot(traj_dict, central_longitude, dynamic_domain):
             dynamic_domain[2] - offset,
             dynamic_domain[3] + offset,
         ],
-        ccrs.PlateCarree(central_longitude=0),
+        #in the following line i substituted 0 with the second central_longitude
+        ccrs.PlateCarree(central_longitude=central_longitude),
     )
 
     # https://stackoverflow.com/questions/65086715/longitude-bounds-of-cartopy-crs-geodetic

src/pytrajplot/plotting/plot_map.py

@@ -15,7 +15,7 @@
 import pandas as pd
 
 # Local
-from ..__init__ import cities_data_path
+from .. import cities_data_path
 from .plot_utils import alps_cities_list
 from .plot_utils import centraleurope_cities_list
 from .plot_utils import ch_cities_list
@@ -216,10 +216,6 @@ def is_visible(lat, lon, domain_boundaries, cross_dateline) -> bool:
                             bool       True if city is within domain boundaries, else false.
 
     """
-    if cross_dateline:
-        if lon < 0:
-            lon = 360 - abs(lon)
-
     lon = float(lon.iloc[0]) if isinstance(lon, pd.Series) else float(lon)
     lat = float(lat.iloc[0]) if isinstance(lat, pd.Series) else float(lat)
     is_in_domain = (
@@ -721,6 +717,64 @@ def get_intersection_point_on_domain_boundaries(
     intersection_point = intersection_point + 0.001 * (point_in - intersection_point)
     return intersection_point
 
+#the next two functions are used later to plot the slice of traj beyond the dateline
+def _unwrap_dateline_series(lon: pd.Series | np.ndarray) -> pd.Series:
+    """Makes the trajectory continuous avoiding the jump -180↔+180 (e.g. -179 → -181)."""
+    arr = np.asarray(lon, dtype=float)
+    if arr.size == 0:
+        return pd.Series(arr, index=getattr(lon, "index", None))
+    diffs = np.diff(arr)
+    offset = 0.0
+    out = [arr[0]]
+    for d, x in zip(diffs, arr[1:]):
+        if d > 180:      # e.g. 179 → -179  (big backwards jump)
+            offset -= 360
+        elif d < -180:   # e.g. -179 → 179  (big forwards jump)
+            offset += 360
+        out.append(x + offset)
+    return pd.Series(out, index=getattr(lon, "index", None))
+
+def create_trajectory_slice_over_dateline(
+    lon: pd.Series | np.ndarray,
+    lat: pd.Series | np.ndarray,
+    threshold: float = 180.0
+) -> tuple[pd.Series, pd.Series]:
+    """
+    gives back the subset (lon_slice, lat_slice) of the points that, after the unwrapping, have |lon| >= threshold (default 180).
+    
+    Parameters
+    ---------
+    lon : pd.Series | np.ndarray
+        Longitudes of the trajectory.
+    lat : pd.Series | np.ndarray
+        Corresponding latitudes (same lenght).
+    threshold : float, default 180.0
+        threshold for the unwrapped longitudes |lon_unwrapped|.
+
+    Returns
+    -------
+    (lon_slice, lat_slice) : tuple[pd.Series, pd.Series]
+        filtered series (maintain the original indexes if lon/lat were Series).
+    """
+    # Convert to Series
+    lon_s = lon if isinstance(lon, pd.Series) else pd.Series(np.asarray(lon, dtype=float))
+    lat_s = lat if isinstance(lat, pd.Series) else pd.Series(np.asarray(lat, dtype=float))
+
+    if len(lon_s) != len(lat_s):
+        raise ValueError("lon e lat must have the same length.")
+
+    # Unwrap longitudes
+    lon_unwrapped = _unwrap_dateline_series(lon_s)
+
+    # identify points over the dateline
+    mask = np.abs(lon_unwrapped.values) >= float(threshold)
+
+    # apply the mask mantaining the index
+    lon_slice = lon_unwrapped[mask]
+    lat_slice = lat_s[mask]
+
+    return lon_slice, lat_slice
+
 
 def add_trajectories_within_domain(
     plot_dict: dict,
@@ -753,10 +807,8 @@ def add_trajectories_within_domain(
             for traj in range(5):
                 latitude = plot_dict["altitude_" + str(i)]["traj_" + str(traj)]["lat"]
                 longitude = plot_dict["altitude_" + str(i)]["traj_" + str(traj)]["lon"]
-                if cross_dateline:
-                    longitude = longitude.apply(
-                        lambda lon: 360 - np.abs(lon) if lon < 0 else lon
-                    )
+                # Unwrap the longitude over the dateline
+                longitude_unwrapped = _unwrap_dateline_series(longitude)
                 ystart = latitude.iloc[0]
                 xstart = longitude.iloc[0]
                 linestyle = subplot_properties_dict[sub_index]
@@ -809,6 +861,21 @@ def add_trajectories_within_domain(
                             transform=ccrs.Geodetic(),
                             rasterized=True,
                         )
+                    # The following plots the slice of the trajectory beyond the dateline that gets clipped
+                    if cross_dateline:
+                        lon_slice, lat_slice= create_trajectory_slice_over_dateline(longitude_unwrapped, latitude, 180.0)
+                        ax.plot(
+                            lon_slice,  # define x-axis
+                            lat_slice,  # define y-axis
+                            linestyle,  # define linestyle
+                            alpha=alpha,  # define line opacity
+                            label=(
+                                None
+                            ),  # We don't want this slice to have its own label
+                            transform=ccrs.Geodetic(),
+                            rasterized=True,
+                        )
+
                     if is_main_trajectory:
                         # add time interval points to main trajectory
                         add_time_interval_points_within_domain(
@@ -894,7 +961,7 @@ def generate_map_plot(
         domain=domain,
         custom_domain_boundaries=trajectory_expansion,
         origin_coordinates=origin_coordinates,
-    )  # sets extent of map
+    )   
     if not domain_boundaries:
         return ax.text(
             0.5,

tests/test_pytrajplot.sh

@@ -28,7 +28,7 @@ pytrajplot $input_dir_tests/test_hres/2_altitudes $output_dir_tests/test_hres/2_
 pytrajplot $input_dir_tests/test_hres/1_altitudes $output_dir_tests/test_hres/1_altitudes --datatype png --domain europe
 
 echo "Test HRES backward trajectory"
-pytrajplot $input_dir_tests/test_hres/backward $output_dir_tests/test_hres/backward --datatype png --domain europe
+pytrajplot $input_dir_tests/test_hres/backward $output_dir_tests/test_hres/backward --datatype png --domain europe --domain dynamic
 
 echo "Test all domains combined w/ HRES model"
 pytrajplot $input_dir_tests/test_hres/4_altitudes $output_dir_tests/test_hres/4_altitudes --datatype png --domain ch --domain alps --domain centraleurope --domain europe --domain dynamic
@@ -37,7 +37,7 @@ echo "Test HRES dateline crossing"
 pytrajplot $input_dir_tests/test_hres/dateline $output_dir_tests/test_hres/dateline --datatype png --domain dynamic
 
 echo "Test HRES w/o side trajectories and the german case"
-pytrajplot $input_dir_tests/test_hres/dateline $output_dir_tests/test_hres/dateline/german --datatype png --domain europe --language de
+pytrajplot $input_dir_tests/test_hres/dateline $output_dir_tests/test_hres/dateline/german --datatype png --domain dynamic --language de
 
 echo "Test HRES Europe with trajectory crossing of zero longitude from east and then leaving domain"
 # (failed in v1.0.0 due to NaNs in the argument of the numpy max function - v1.0.1 uses nanmax instead)

tests/test_recent_opr.sh

@@ -6,14 +6,14 @@
 # Conda default environment
 conda_env=pytrajplot
 # Lagranto-output storage location
-store_osm=/store/mch/msopr/osm
+store_osm=/store_new/mch/msopr/osm
 # pytrajplot output location
-pytrajplot_out=.
+pytrajplot_out=local
 # Graphics format
-datatype_opt="--datatype png --datatype pdf"
+datatype_opt="--datatype png" # --datatype pdf"
 # Domain options for
-# COSMO-1E-CTR (c1), IFS-HRES-Europe (ie), IFS-HRES global (ig)
-domain_opt_c1="--domain ch --domain alps"
+# COSMO-1E-CTR (i1), IFS-HRES-Europe (ie), IFS-HRES global (ig)
+domain_opt_i1="--domain ch --domain alps"
 domain_opt_ie="--domain alps --domain centraleurope --domain europe"
 domain_opt_ig="--domain dynamic --domain dynamic_zoom"
 # --------
@@ -32,7 +32,8 @@ bt_00=$(date --utc --date="today 00" +%Y%m%d%H)
 bt_03=$(date --utc --date="today 03" +%Y%m%d%H)
 
 # Load conda env for pytrajplot if CONDA_PREFIX not defined
-[[ -z $CONDA_PREFIX ]] && conda activate $conda_env
+#[[ -z $CONDA_PREFIX ]] &&
+conda activate $conda_env
 echo CONDA_PREFIX=$CONDA_PREFIX
 # Report version
 $CONDA_PREFIX/bin/pytrajplot --version
@@ -46,17 +47,17 @@ for basetime in $bt_06 $bt_12 $bt_18 $bt_21 $bt_00 $bt_03 ; do
     echo Basetime: $(date --utc --date="${basetime:0:8} $hh" "+%F %H UTC")
 
     # Operational INPUT_DIRs:
-    input_dir_c1=$(echo $store_osm/COSMO-1E/FCST${yy}/${yymmddhh}_4??/lagranto_c/000)
+    input_dir_i1=$(echo $store_osm/ICON-CH1-EPS/FCST${yy}/${yymmddhh}_???/lagranto_c/000)
     input_dir_ie=$store_osm/IFS-HRES/IFS-HRES-LAGRANTO${yy}/${yymmddhh}_LIH/lagranto_f
     input_dir_ig=$store_osm/IFS-HRES/IFS-HRES-LAGRANTO${yy}/${yymmddhh}_LIH/lagranto_c
 
     # Output directories
-    output_dir_c1=$pytrajplot_out/plot_c1_${basetime}
+    output_dir_i1=$pytrajplot_out/plot_i1_${basetime}
     output_dir_ie=$pytrajplot_out/plot_ie_${basetime}
     output_dir_ig=$pytrajplot_out/plot_ig_${basetime}
 
     # Submit jobs
-    for model in c1 ie ig ; do
+    for model in i1 ie ig ; do
         eval input_dir=\$input_dir_$model
         eval output_dir=\$output_dir_$model
         eval domain_opt=\$domain_opt_$model