From 5642aa52a10a892f4362bdf36f506a61d24cacb9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Peider=20K=C3=B6nz?= Date: Mon, 6 Oct 2025 10:54:21 +0200 Subject: [PATCH 1/5] general corrections of the code --- README.md | 7 +++---- src/pytrajplot/parse_data.py | 24 ++++++++++++------------ tests/test_pytrajplot.sh | 4 ++-- tests/test_recent_opr.sh | 19 ++++++++++--------- 4 files changed, 27 insertions(+), 27 deletions(-) diff --git a/README.md b/README.md index 35e4e08..41f5ba4 100644 --- a/README.md +++ b/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//bin/pip`). diff --git a/src/pytrajplot/parse_data.py b/src/pytrajplot/parse_data.py index b5927c1..4649856 100644 --- a/src/pytrajplot/parse_data.py +++ b/src/pytrajplot/parse_data.py @@ -360,8 +360,8 @@ 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 @@ -369,25 +369,25 @@ def read_trajectory(trajectory_file_path, start_df, plot_info_dict): # 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 diff --git a/tests/test_pytrajplot.sh b/tests/test_pytrajplot.sh index 32bec31..40f126a 100755 --- a/tests/test_pytrajplot.sh +++ b/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) diff --git a/tests/test_recent_opr.sh b/tests/test_recent_opr.sh index e197c83..c36264b 100755 --- a/tests/test_recent_opr.sh +++ b/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 From 51d7f314411db1e5d85a9c27e834faebf63b099d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Peider=20K=C3=B6nz?= Date: Mon, 6 Oct 2025 10:57:44 +0200 Subject: [PATCH 2/5] fix bug concerning trajectories crossing the dateline --- .../plotting/analyse_trajectories.py | 5 +- src/pytrajplot/plotting/plot_map.py | 87 ++++++++++++++++--- 2 files changed, 80 insertions(+), 12 deletions(-) diff --git a/src/pytrajplot/plotting/analyse_trajectories.py b/src/pytrajplot/plotting/analyse_trajectories.py index 46cb8b3..3314e75 100644 --- a/src/pytrajplot/plotting/analyse_trajectories.py +++ b/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 diff --git a/src/pytrajplot/plotting/plot_map.py b/src/pytrajplot/plotting/plot_map.py index 7d301c7..5b34385 100644 --- a/src/pytrajplot/plotting/plot_map.py +++ b/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 sre 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 - ) + #NEW: unwraps 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, ) + #I ADDED THE FOLLOWING IF: IT PLOTS THE SLICE OF TRAJ CLIPPED OVER THE DATELINE + 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, From 35a88644bc49b53ae3eb7f18d7560719725cbff5 Mon Sep 17 00:00:00 2001 From: Peiderk Date: Thu, 9 Oct 2025 11:44:45 +0200 Subject: [PATCH 3/5] Update src/pytrajplot/plotting/plot_map.py Co-authored-by: Pirmin Kaufmann --- src/pytrajplot/plotting/plot_map.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pytrajplot/plotting/plot_map.py b/src/pytrajplot/plotting/plot_map.py index 5b34385..483278e 100644 --- a/src/pytrajplot/plotting/plot_map.py +++ b/src/pytrajplot/plotting/plot_map.py @@ -717,7 +717,7 @@ def get_intersection_point_on_domain_boundaries( intersection_point = intersection_point + 0.001 * (point_in - intersection_point) return intersection_point -#the next two functions sre used later to plot the slice of traj beyond the dateline +#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) From 35228e4b9ef6c06e7aa0c946006bf70477b5bbbd Mon Sep 17 00:00:00 2001 From: Peiderk Date: Thu, 9 Oct 2025 11:44:52 +0200 Subject: [PATCH 4/5] Update src/pytrajplot/plotting/plot_map.py Co-authored-by: Pirmin Kaufmann --- src/pytrajplot/plotting/plot_map.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pytrajplot/plotting/plot_map.py b/src/pytrajplot/plotting/plot_map.py index 483278e..eb3bcef 100644 --- a/src/pytrajplot/plotting/plot_map.py +++ b/src/pytrajplot/plotting/plot_map.py @@ -861,7 +861,7 @@ def add_trajectories_within_domain( transform=ccrs.Geodetic(), rasterized=True, ) - #I ADDED THE FOLLOWING IF: IT PLOTS THE SLICE OF TRAJ CLIPPED OVER THE DATELINE + # 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( From 9d7815e5b511a6e75cffcb4a20674d2691b3373f Mon Sep 17 00:00:00 2001 From: Peiderk Date: Thu, 9 Oct 2025 11:44:59 +0200 Subject: [PATCH 5/5] Update src/pytrajplot/plotting/plot_map.py Co-authored-by: Pirmin Kaufmann --- src/pytrajplot/plotting/plot_map.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pytrajplot/plotting/plot_map.py b/src/pytrajplot/plotting/plot_map.py index eb3bcef..2fd876a 100644 --- a/src/pytrajplot/plotting/plot_map.py +++ b/src/pytrajplot/plotting/plot_map.py @@ -807,7 +807,7 @@ 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"] - #NEW: unwraps the longitude over the dateline + # Unwrap the longitude over the dateline longitude_unwrapped = _unwrap_dateline_series(longitude) ystart = latitude.iloc[0] xstart = longitude.iloc[0]