Add option to preserve 0–360 longitude convention (--keep_longitude)

inputs/namelist

@@ -1,12 +1,11 @@
 ;standard_name;Variable;Units
-Air Temperature;air_temperature;TMP_2_ISBL;K
-Geopotential Height;geopotential_height;HGT_2_ISBL;m
-Specific Humidity;specific_humidity;q_2_ISBL;kg/kg
-Omega Velocity;omega;V_VEL_2_ISBL;Pa/s
-Eastward Wind Component;eastward_wind;U_GRD_2_ISBL;m/s
-Northward Wind Component;northward_wind;V_GRD_2_ISBL;m/s
-Longitude;;lon_2
-Latitude;;lat_2
-Time;;initial_time0_hours
-Vertical Level;;lv_ISBL3
-
+Air Temperature;air_temperature;t;K
+Geopotential;geopotential;z;m**2/s**2
+Specific Humidity;specific_humidity;q;kg/kg
+Omega Velocity;omega;w;Pa/s
+Eastward Wind Component;eastward_wind;u;m/s
+Northward Wind Component;northward_wind;v;m/s
+Longitude;;longitude
+Latitude;;latitude
+Time;;valid_time
+Vertical Level;;pressure_level

inputs/namelist_ERA5

@@ -1,11 +1,11 @@
 ;standard_name;Variable;Units
-Air Temperature;air_temperature;T;K
-Geopotential;geopotential;Z;m**2/s**2
-Specific Humidity;specific_humidity;Q;kg/kg
-Omega Velocity;omega;W;Pa/s
-Eastward Wind Component;eastward_wind;U;m/s
-Northward Wind Component;northward_wind;V;m/s
+Air Temperature;air_temperature;t;K
+Geopotential;geopotential;z;m**2/s**2
+Specific Humidity;specific_humidity;q;kg/kg
+Omega Velocity;omega;w;Pa/s
+Eastward Wind Component;eastward_wind;u;m/s
+Northward Wind Component;northward_wind;v;m/s
 Longitude;;longitude
 Latitude;;latitude
-Time;;time
-Vertical Level;;level
+Time;;valid_time
+Vertical Level;;pressure_level

src/cli_interface.py

@@ -35,6 +35,10 @@ def parse_arguments(custom_args=None):
     parser.add_argument("infile", 
                         help="""Input .nc file with temperature, specific humidity and meridional, zonal, and vertical components of the wind, in pressure levels.""")
 
+    parser.add_argument("--keep_longitude",
+                        action="store_true",
+                        help="Keep original longitude convention from input file. If not used, longitudes are converted from 0–360 to -180–180.")
+
     group = parser.add_mutually_exclusive_group(required=True)
     group.add_argument("-f", "--fixed", action='store_true',
                        help="Compute the energetics for a Fixed domain specified by the 'inputs/box_limits' file.")
@@ -63,4 +67,4 @@ def parse_arguments(custom_args=None):
     else:
         args = parser.parse_args()
 
-    return args
+    return args

src/data_handling.py

@@ -37,19 +37,21 @@ def load_data(infile, longitude_indexer, args, app_logger):
         app_logger.info(f'⏳ Loading {infile}...')
         with dask.config.set(array={'slicing': {'split_large_chunks': True}}):
             if args.gfs:
-                data = convert_lon(
-                    xr.open_mfdataset(
+                data = xr.open_mfdataset(
                         infile,
                         engine='cfgrib',
                         parallel=True,
                         filter_by_keys={'typeOfLevel': 'isobaricInhPa'},
                         combine='nested',
                         concat_dim='time'
-                    ),
-                    longitude_indexer
-                )
+                    )
             else:
-                data = convert_lon(xr.open_dataset(infile), longitude_indexer)
+                data = xr.open_dataset(infile)
+
+            if not getattr(args, "keep_longitude", False):
+                data = convert_lon(data, longitude_indexer)
+            else:
+                app_logger.info('Keeping original longitude convention from input file')
 
         app_logger.info(f'✅ Loaded {infile} successfully!')
         return data

src/select_domain.py

@@ -31,7 +31,24 @@
 import matplotlib.ticker as ticker
 
 nclicks = 2
-crs_longlat = ccrs.PlateCarree() 
+
+# CRS used by the input data.
+# This should remain PlateCarree() because the data coordinates themselves
+# are still regular latitude/longitude coordinates.
+data_crs = ccrs.PlateCarree()
+
+def get_map_crs(args=None):
+    """
+    Return the map projection.
+
+    If --keep_longitude is used, the input data may remain in 0–360 longitude.
+    In this case, centering the map at 180 degrees avoids dateline-related
+    plotting artifacts in Cartopy.
+    """
+    if args is not None and getattr(args, "keep_longitude", False):
+        return ccrs.PlateCarree(central_longitude=180)
+
+    return ccrs.PlateCarree()
 
 # Transformation function
 def coordXform(orig_crs, target_crs, x, y):
@@ -73,7 +90,7 @@ def fmt(x, pos):
     b = int(b)
     return r'${} \times 10^{{{}}}$'.format(a, b)
 
-def draw_box(ax, limits, crs):
+def draw_box(ax, limits, crs=data_crs):
     """
     Draws a rectangular box on the plot based on specified limits.
 
@@ -87,7 +104,7 @@ def draw_box(ax, limits, crs):
     pgon = Polygon([(min_lon, min_lat), (min_lon, max_lat), (max_lon, max_lat), (max_lon, min_lat), (min_lon, min_lat)])
     ax.add_geometries([pgon], crs=crs, facecolor='None', edgecolor='k', linewidth=3, alpha=1, zorder=3)
 
-def plot_zeta(ax, zeta, lat, lon, hgt=None):
+def plot_zeta(ax, zeta, lat, lon, hgt=None, crs=data_crs):
     """
     Plots the vorticity field and optionally the geopotential height.
 
@@ -100,10 +117,10 @@ def plot_zeta(ax, zeta, lat, lon, hgt=None):
     """
     norm = colors.TwoSlopeNorm(vmin=min(-zeta.max(), zeta.min()), vcenter=0, vmax=max(zeta.max(), -zeta.min()))
     cmap = cmo.balance
-    cf1 = ax.contourf(lon, lat, zeta, cmap=cmap, norm=norm, levels=51, transform=crs_longlat) 
+    cf1 = ax.contourf(lon, lat, zeta, cmap=cmap, norm=norm, levels=51, transform=data_crs) 
     plt.colorbar(cf1, pad=0.1, orientation='vertical', shrink=0.5, format=ticker.FuncFormatter(fmt))
     if hgt is not None:
-        cs = ax.contour(lon, lat, hgt, levels=11, colors='#747578', linestyles='dashed', linewidths=1, transform=crs_longlat)
+        cs = ax.contour(lon, lat, hgt, levels=11, colors='#747578', linestyles='dashed', linewidths=1, transform=crs)
         ax.clabel(cs, cs.levels, inline=True, fontsize=10)
 
 def map_decorators(ax):
@@ -150,18 +167,18 @@ def plot_min_max_zeta(ax, zeta, lat, lon, limits, args):
                 for point in points:
                     ax.scatter(point[lon.dims[0]], point[lat.dims[0]],
                            marker='o', facecolors='none', linewidth=3,
-                           edgecolor='k',  s=200)
+                           edgecolor='k',  s=200, transform=data_crs)
         else:
             for point in min_max_zeta_loc:
                 ax.scatter(point[lon.dims[0]], point[lat.dims[0]],
                        marker='o', facecolors='none', linewidth=3,
-                       edgecolor='k',  s=200)
+                       edgecolor='k',  s=200, transform=data_crs)
     else:
         ax.scatter(min_max_zeta_loc[lon.dims[0]], min_max_zeta_loc[lat.dims[0]],
                marker='o', facecolors='none', linewidth=3,
-               edgecolor='k',  s=200)
+               edgecolor='k',  s=200, transform=data_crs)
 
-def initial_domain(zeta, lat, lon):
+def initial_domain(zeta, lat, lon, args=None):
     """
     Interactively selects an initial spatial domain on a map.
 
@@ -175,10 +192,11 @@ def initial_domain(zeta, lat, lon):
     """
     plt.close('all')
     fig = plt.figure(figsize=(10, 8))
-    ax = plt.axes(projection=crs_longlat)
+    map_crs = get_map_crs(args)
+    ax = plt.axes(projection=map_crs)
     fig.add_axes(ax)
     ax.set_global()
-    plot_zeta(ax, zeta, lat, lon)
+    plot_zeta(ax, zeta, lat, lon, crs=data_crs)
     map_decorators(ax)
     plt.subplots_adjust(bottom=0, top=1.2)
 
@@ -195,13 +213,15 @@ def initial_domain(zeta, lat, lon):
         xmin,xmax = min([pts[0,0],pts[1,0]]),max([pts[0,0],pts[1,0]])
         ymin,ymax = min([pts[0,1],pts[1,1]]),max([pts[0,1],pts[1,1]])
         xs, ys = np.array((xmin,xmax)), np.array((ymin,ymax))
-        lls = coordXform(crs_longlat, crs_longlat, xs, ys)
+        lls = coordXform(map_crs, data_crs, xs, ys)
+        if args is not None and getattr(args, "keep_longitude", False):
+            lls[:, 0] = lls[:, 0] % 360
         min_lon,min_lat = lls[0,0], lls[0,1]
         max_lon, max_lat = lls[1,0], lls[1,1]
 
         limits = {'min_lon':min_lon,'max_lon':max_lon,
                 'min_lat':min_lat, 'max_lat':max_lat}
-        draw_box(ax, limits, crs_longlat)
+        draw_box(ax, limits, data_crs)
 
         tellme('Happy? Key press any keyboard key for yes, mouse click for no')
         if plt.waitforbuttonpress():
@@ -228,12 +248,13 @@ def draw_box_map(u, v, zeta, hgt, lat, lon, timestr, args):
     """
     plt.close('all')
     fig = plt.figure(figsize=(10, 8))
-    ax = plt.axes(projection=crs_longlat)
+    map_crs = get_map_crs(args)
+    ax = plt.axes(projection=map_crs)
     fig.add_axes(ax)
 
-    plot_zeta(ax, zeta, lat, lon, hgt)
+    plot_zeta(ax, zeta, lat, lon, hgt, crs=data_crs)
     ax.streamplot(lon.values, lat.values, u.values, v.values, color='#2A1D21',
-              transform=crs_longlat)
+              transform=data_crs)
     map_decorators(ax)
 
     # Convert to datetime object
@@ -257,13 +278,15 @@ def draw_box_map(u, v, zeta, hgt, lat, lon, timestr, args):
         xmin,xmax = min([pts[0,0],pts[1,0]]),max([pts[0,0],pts[1,0]])
         ymin,ymax = min([pts[0,1],pts[1,1]]),max([pts[0,1],pts[1,1]])
         xs, ys = np.array((xmin,xmax)), np.array((ymin,ymax))
-        lls = coordXform(crs_longlat, crs_longlat, xs, ys)
+        lls = coordXform(map_crs, data_crs, xs, ys)
+        if getattr(args, "keep_longitude", False):
+            lls[:, 0] = lls[:, 0] % 360
         min_lon,min_lat = lls[0,0], lls[0,1]
         max_lon, max_lat = lls[1,0], lls[1,1]
 
         limits = {'min_lon':min_lon,'max_lon':max_lon,
                 'min_lat':min_lat, 'max_lat':max_lat}
-        draw_box(ax, limits, crs_longlat)
+        draw_box(ax, limits, data_crs)
         plot_min_max_zeta(ax, zeta, lat, lon, limits, args)
 
         tellme('Happy? Key press any keyboard key for yes, mouse click for no')
@@ -344,4 +367,4 @@ def get_domain_limits(args, *variables_at_plevel, track=None):
             'central_lon': (max_lon + min_lon) / 2
             }
 
-    return current_domain_limits
+    return current_domain_limits

src/utils.py

@@ -61,6 +61,7 @@ def initialize_logging(results_subdirectory, args):
 
     return app_logger
 
+
 def convert_lon(input_data, longitude_indexer):
     """
 
@@ -83,6 +84,7 @@ def convert_lon(input_data, longitude_indexer):
     input_data = input_data.sortby(input_data[longitude_indexer])
     return input_data
 
+
 def handle_track_file(input_data, times, longitude_indexer, latitude_indexer, app_logger):
     """
     Handles the track file by validating its time and spatial limits against the provided dataset.
@@ -270,7 +272,7 @@ def slice_domain(input_data, args, namelist_df):
         zeta = vorticity(iu_plevel, iv_plevel).metpy.dequantify()
 
         lat, lon = iu_plevel[latitude_indexer], iu_plevel[longitude_indexer]
-        domain_limits = initial_domain(zeta, lat, lon)
+        domain_limits = initial_domain(zeta, lat, lon, args)
         WesternLimit = domain_limits['min_lon']
         EasternLimit = domain_limits['max_lon']
         SouthernLimit = domain_limits['min_lat']
@@ -333,4 +335,4 @@ def get_domain_extreme_values(itime, args, slices_plevel, track=None):
         min_max_hgt = float(ight_plevel_slice.min())
         max_wind = float(iwspd_plevel_slice.max())
 
-    return min_max_zeta, min_max_hgt, max_wind
+    return min_max_zeta, min_max_hgt, max_wind

src/visualization.py

@@ -41,6 +41,25 @@
 from cartopy.feature import BORDERS
 from src.select_domain import plot_zeta 
 
+def get_cartopy_crs(args):
+    """
+    Return data and map coordinate reference systems.
+
+    data_crs is always PlateCarree because the data are stored as regular
+    latitude/longitude coordinates.
+
+    map_crs changes only when --keep_longitude is used. In that case, the map
+    is centered at 180° to correctly display data kept in 0–360 longitude.
+    """
+    data_crs = ccrs.PlateCarree()
+
+    if getattr(args, "keep_longitude", False):
+        map_crs = ccrs.PlateCarree(central_longitude=180)
+    else:
+        map_crs = ccrs.PlateCarree()
+
+    return data_crs, map_crs
+
 def map_features(ax):
     """
     Adds coastline and border features to the matplotlib Axes.
@@ -95,17 +114,19 @@ def plot_fixed_domain(limits, data_plevel, args, results_subdirectory, time, app
 
     # Create figure
     plt.close('all')
-    fig, ax = plt.subplots(figsize=(8, 8.5), subplot_kw=dict(projection=ccrs.PlateCarree()))
+
+    data_crs, map_crs = get_cartopy_crs(args)
+    fig, ax = plt.subplots(figsize=(8, 8.5), subplot_kw=dict(projection=map_crs))
 
     # Set map extent and features
-    ax.set_extent([min_lon-20, max_lon+20, max_lat+20, min_lat-20], crs=ccrs.PlateCarree())
+    ax.set_extent([min_lon-20, max_lon+20, max_lat+20, min_lat-20], crs=data_crs)
     map_features(ax)
     Brazil_states(ax, facecolor='None')
 
     # Plot selected domain
     # Create a sample polygon, `pgon`
     pgon = Polygon(((min_lon, min_lat), (min_lon, max_lat), (max_lon, max_lat), (max_lon, min_lat), (min_lon, min_lat)))
-    ax.add_geometries([pgon], crs=ccrs.PlateCarree(), facecolor='None', edgecolor='#BF3D3B', linewidth=3, alpha=1, zorder=3)
+    ax.add_geometries([pgon], crs=data_crs, facecolor='None', edgecolor='#BF3D3B', linewidth=3, alpha=1, zorder=3)
 
     # Add gridlines
     gl = ax.gridlines(draw_labels=True,zorder=2)    
@@ -123,13 +144,13 @@ def plot_fixed_domain(limits, data_plevel, args, results_subdirectory, time, app
     Brazil_states(ax, facecolor='None')
 
     # Plot central point, mininum vorticity, minimum hgt and maximum wind
-    ax.scatter(central_lon, central_lat,  marker='o', c='#31332e', s=100, zorder=4)
+    ax.scatter(central_lon, central_lat,  marker='o', c='#31332e', s=100, zorder=4, transform=data_crs)
     ax.scatter(data_plevel[f'{args.track_vorticity}_zeta_{args.level}']['longitude'], data_plevel[f'{args.track_vorticity}_zeta_{args.level}']['latitude'],
-                marker='s', c='#31332e', s=100, zorder=4, label=f'{args.track_vorticity} zeta')
+                marker='s', c='#31332e', s=100, zorder=4, label=f'{args.track_vorticity} zeta', transform=data_crs)
     ax.scatter(data_plevel[f'{args.track_geopotential}_hgt_{args.level}']['longitude'], data_plevel[f'{args.track_geopotential}_hgt_{args.level}']['latitude'],
-                marker='x', c='#31332e', s=100, zorder=4, label=f'{args.track_geopotential} hgt')
+                marker='x', c='#31332e', s=100, zorder=4, label=f'{args.track_geopotential} hgt', transform=data_crs)
     ax.scatter(data_plevel['max_wind']['longitude'], data_plevel['max_wind']['latitude'],
-                marker='^', c='#31332e', s=100, zorder=4, label='max wind')
+                marker='^', c='#31332e', s=100, zorder=4, label='max wind', transform=data_crs)
     plt.legend(loc='upper left', frameon=True, fontsize=14, bbox_to_anchor=(1.1,1.2))
 
     # Save figure
@@ -167,10 +188,12 @@ def plot_track(track, args, figures_directory, app_logger):
         min_lat, max_lat = track['Lat'].min(), track['Lat'].max()
 
         plt.close('all')
+
+        data_crs, map_crs = get_cartopy_crs(args)
         fig, ax = plt.subplots(figsize=(12, 9) if (max_lon - min_lon) <= (max_lat - min_lat) else (14, 9), 
-                               subplot_kw={"projection": ccrs.PlateCarree()})
+                               subplot_kw={"projection": map_crs})
 
-        ax.set_extent([min_lon - 10, max_lon + 10, min_lat - 10, max_lat + 10], crs=ccrs.PlateCarree())
+        ax.set_extent([min_lon - 10, max_lon + 10, min_lat - 10, max_lat + 10], crs=data_crs)
         ax.coastlines()
         ax.add_feature(cartopy.feature.LAND)
         ax.add_feature(cartopy.feature.OCEAN, facecolor="lightblue")
@@ -179,19 +202,19 @@ def plot_track(track, args, figures_directory, app_logger):
         gl.top_labels = gl.right_labels = False
 
         # Plotting the track
-        ax.plot(track['Lon'], track['Lat'], color='#383838')
+        ax.plot(track['Lon'], track['Lat'], color='#383838', transform=data_crs)
 
         # Enhanced visualization with normalized wind speed and vorticity
         if f'{args.track_vorticity}_zeta_{args.level}' in track and f'max_wind_{args.level}' in track:
             normalized_sizes = normalize(track[[f'max_wind_{args.level}']].values.reshape(1, -1))**4 * 100000
             scatter = ax.scatter(track['Lon'], track['Lat'], zorder=3, c=track[f'{args.track_vorticity}_zeta_{args.level}'].values, 
                                  cmap=cmo.deep_r, edgecolor='gray', s=normalized_sizes.flatten(), 
-                                 label=f'Normalized max wind speed at {args.level} hPa')
+                                 label=f'Normalized max wind speed at {args.level} hPa', transform=data_crs)
             plt.colorbar(scatter, pad=0.1, orientation='vertical', shrink=0.5, label=f'{args.track_vorticity} vorticity')
 
         # Marking the start and end points
-        ax.text(track.iloc[0]['Lon'], track.iloc[0]['Lat'], 'A', fontsize=12, ha='center', va='center', color='green')
-        ax.text(track.iloc[-1]['Lon'], track.iloc[-1]['Lat'], 'Z', fontsize=12, ha='center', va='center', color='red')
+        ax.text(track.iloc[0]['Lon'], track.iloc[0]['Lat'], 'A', fontsize=12, ha='center', va='center', color='green', transform=data_crs)
+        ax.text(track.iloc[-1]['Lon'], track.iloc[-1]['Lat'], 'Z', fontsize=12, ha='center', va='center', color='red', transform=data_crs)
 
         filename = os.path.join(figures_directory, 'track_boxes.png')
         plt.savefig(filename, bbox_inches='tight')
@@ -303,4 +326,4 @@ def hovmoller_mean_zeta(Zeta, figures_subdirectory, app_logger):
         app_logger.info(f'🌍 Hovmöller diagram of mean zeta created and saved: {filename}')
 
     except Exception as e:
-        app_logger.error(f"❌ Failed to create Hovmöller diagram of mean zeta: {e}")
+        app_logger.error(f"❌ Failed to create Hovmöller diagram of mean zeta: {e}")