diff --git a/inputs/Extract_coords_time-lat-lon-length-width.py b/inputs/Extract_coords_time-lat-lon-length-width.py new file mode 100644 index 0000000..866fc76 --- /dev/null +++ b/inputs/Extract_coords_time-lat-lon-length-width.py @@ -0,0 +1,399 @@ +#!/usr/bin/env python3 +""" +Interactive cyclone tracking tool. + +This script reads a NetCDF file, displays a selected meteorological field +for each time step, and allows the user to draw a rectangular box around +a cyclone using the mouse. + +For each selected box, the script saves: + - time + - centroid latitude + - centroid longitude + - box length in degrees longitude + - box width in degrees latitude + +The output file follows the format: + + time;Lat;Lon;length;width + +Example: + 2026-04-01-0000;-35.2500;190.5000;8.0000;6.5000 + +The script preserves the longitude convention of the input NetCDF file: + - 0 to 360 + - -180 to 180 + +Optionally, a continent shapefile can be overlaid using pyshp. +""" + +import os +import argparse + +import xarray as xr +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.widgets import RectangleSelector +import shapefile + +# ========================================================== +# Configurações +# ========================================================== +parser = argparse.ArgumentParser( + description=( + "Interactively extract cyclone centroid coordinates and box " + "dimensions from a NetCDF file." + ) +) + +parser.add_argument( + "file_nc", + type=str, + help="Path to the input NetCDF file." +) + +parser.add_argument( + "--var", + type=str, + default="msl", + help="Variable name in the NetCDF file. Default: msl" +) + +parser.add_argument( + "--level", + type=float, + default=None, + help=( + "Optional pressure level to select if the dataset has a vertical " + "pressure-level coordinate. Example: --level 850." + ) +) + +parser.add_argument( + "--out", + type=str, + default="track_ciclone.txt", + help="Name of the output file (default: track_ciclone.txt)" +) + +parser.add_argument("--shp", type=str, default=None, + help="Optional continent shapefile path.") + +args = parser.parse_args() + +file_nc = args.file_nc +outfile = args.out +varname = args.var +level_value = args.level +shp_path = args.shp + +# ========================================================== +# Input checks +# ========================================================== +if not os.path.exists(file_nc): + raise FileNotFoundError(f"File not found: {file_nc}") + +if shp_path and not os.path.exists(shp_path): + raise FileNotFoundError(f"Shapefile not found: {shp_path}") + +# ========================================================== +# Open dataset +# ========================================================== +ds = xr.open_dataset(file_nc) + +# ========================================================== +# Select pressure level, if requested +# ========================================================== +if level_value is not None: + possible_level_names = [ + "level", + "pressure_level", + "isobaricInhPa", + "plev", + "lev" + ] + + level_name = None + + for name in possible_level_names: + if name in ds.coords or name in ds.dims: + level_name = name + break + + if level_name is None: + raise ValueError( + "The --level option was used, but no pressure-level coordinate " + "was found in the NetCDF file." + ) + + ds = ds.sel({level_name: level_value}, method="nearest") + + print(f"Selected level: {level_name} = {float(ds[level_name].values)}") + +# ========================================================== +# Detect coordinate names automatically +# ========================================================== +lat_name = None +lon_name = None +time_name = None + +for name in ds.coords: + lname = name.lower() + if lname in ["lat", "latitude"]: + lat_name = name + elif lname in ["lon", "longitude"]: + lon_name = name + elif lname in ["time", "valid_time", "date"]: + time_name = name + +if lat_name is None or lon_name is None or time_name is None: + raise ValueError( + "Could not automatically detect latitude, longitude, or time " + "coordinates in the NetCDF file." + ) + +# ========================================================== +# Domain and longitude convention +# ========================================================== +lat_min = float(ds[lat_name].min()) +lat_max = float(ds[lat_name].max()) +lon_min = float(ds[lon_name].min()) +lon_max = float(ds[lon_name].max()) + +use_lon_360 = lon_max > 180 + +print( + "Detected longitude convention: " + f"{'0-360' if use_lon_360 else '-180 to 180'}" +) + +# ========================================================== +# Shapefile plotting function +# ========================================================== +def plot_shapefile(ax, + shp_path, + use_lon_360, + edgecolor="black", + linewidth=0.6, + label=None): + """ + Plot shapefile contours using Matplotlib. + + Parameters + ---------- + ax : matplotlib.axes.Axes + Axis where the shapefile contours will be drawn. + shp_path : str + Path to the shapefile. + use_lon_360 : bool + If True, shapefile longitudes are converted from -180/180 to 0/360. + If False, shapefile longitudes are kept or converted to -180/180. + edgecolor : str + Contour color. + linewidth : float + Contour line width. + label : str or None + Optional label for the legend. The label is added only once. + + Notes + ----- + The function splits line segments when large longitude jumps are detected. + This avoids artificial lines near the dateline. + """ + + sf = shapefile.Reader(shp_path) + first = True + + for shp in sf.shapes(): + pts = np.array(shp.points) + if pts.size == 0: + continue + + x = pts[:, 0] + y = pts[:, 1] + + if use_lon_360: + x = np.where(x < 0, x + 360, x) + else: + x = np.where(x > 180, x - 360, x) + + # Split multipart geometries. + parts = list(shp.parts) + [len(pts)] + + for i in range(len(parts) - 1): + xs = x[parts[i]:parts[i+1]] + ys = y[parts[i]:parts[i+1]] + + # Split segments that cross the dateline after longitude conversion. + jumps = np.abs(np.diff(xs)) > 180 + + start = 0 + + for j, jump in enumerate(jumps): + if jump: + if first and label is not None: + ax.plot(xs[start:j+1], ys[start:j+1], + color=edgecolor, linewidth=linewidth, label=label) + first = False + else: + ax.plot(xs[start:j+1], ys[start:j+1], + color=edgecolor, linewidth=linewidth) + start = j + 1 + + if first and label is not None: + ax.plot(xs[start:], ys[start:], color=edgecolor, linewidth=linewidth, label=label) + first = False + else: + ax.plot(xs[start:], ys[start:], color=edgecolor, linewidth=linewidth) + +# ========================================================== +# Prepare output file +# time;Lat;Lon,length;width +# ========================================================== +nt = len(ds[time_name]) +da_all = ds[varname] + +with open(outfile, "w") as f: + f.write("time;Lat;Lon;length;width\n") + +# ========================================================== +# Time loop +# ========================================================== +for n in range(nt): + + da = da_all.isel({time_name: n}) + + # Convert Pa to hPa when the field appears to be pressure in Pa. + if float(da.mean()) > 2000: + da_plot = da / 100.0 + else: + da_plot = da + + time_value = pd.to_datetime(ds[time_name].isel({time_name: n}).values) + time_txt = time_value.strftime("%Y-%m-%d-%H%M") + time_title = time_value.strftime("%HZ%d%b%Y").upper() + + # ====================================================== + # Create plot + # ====================================================== + fig, ax = plt.subplots(figsize=(7, 9)) + + vmin = np.floor(float(da_plot.min())) + vmax = np.ceil(float(da_plot.max())) + levels = np.arange(vmin, vmax + 4, 4) + + cs = ax.contour( + ds[lon_name].values, + ds[lat_name].values, + da_plot.values, + levels=levels, + colors="black", + linewidths=1 + ) + + ax.clabel(cs, inline=True, fontsize=8, fmt="%.0f") + + ax.set_xlim(lon_min, lon_max) + ax.set_ylim(lat_min, lat_max) + ax.set_xlabel("Longitude") + ax.set_ylabel("Latitude") + ax.set_title(f"Mean sea level pressure {time_title}", fontsize=14) + + ax.grid(True, linestyle="--", alpha=0.4) + + # Plot optional continent shapefile + if shp_path is not None: + plot_shapefile( + ax, + shp_path, + use_lon_360, + edgecolor="red", + linewidth=1.5, + label="Continents" + ) + ax.legend(loc="upper right") + + print("\n======================================") + print(f"Time step {n+1}/{nt}: {time_txt}") + print("Left-click, drag the rectangle, and release to select the cyclone box.") + print("Close the window to skip this time step.") + print("======================================") + + selection = {} + + def onselect(eclick, erelease): + """ + Callback executed after drawing the rectangle. + + The selected box limits are stored in the 'selection' dictionary. + """ + lon_a, lat_a = eclick.xdata, eclick.ydata + lon_b, lat_b = erelease.xdata, erelease.ydata + + if lon_a is None or lon_b is None or lat_a is None or lat_b is None: + return + + lon_left = min(lon_a, lon_b) + lon_right = max(lon_a, lon_b) + lat_bottom = min(lat_a, lat_b) + lat_top = max(lat_a, lat_b) + + selection["lon_left"] = lon_left + selection["lon_right"] = lon_right + selection["lat_bottom"] = lat_bottom + selection["lat_top"] = lat_top + + plt.close(fig) + + selector = RectangleSelector( + ax, + onselect, + useblit=True, + button=[1], + minspanx=0.1, + minspany=0.1, + spancoords="data", + interactive=True + ) + + plt.show() + + if not selection: + print("No box selected. Skipping this time step.") + plt.close(fig) + continue + + lon_left = selection["lon_left"] + lon_right = selection["lon_right"] + lat_bottom = selection["lat_bottom"] + lat_top = selection["lat_top"] + + # Compute box centroid and dimensions in degrees. + lon_center = (lon_left + lon_right) / 2.0 + lat_center = (lat_bottom + lat_top) / 2.0 + + length = lon_right - lon_left + width = lat_top - lat_bottom + + # Save selected box information. + with open(outfile, "a") as f: + f.write( + f"{time_txt};" + f"{lat_center:.4f};" + f"{lon_center:.4f};" + f"{length:.4f};" + f"{width:.4f}\n" + ) + + print( + f"Saved: {time_txt};" + f"{lat_center:.4f};" + f"{lon_center:.4f};" + f"{length:.4f};" + f"{width:.4f}" + ) + + +print("\nProcessing completed.") +print(f"Output file saved as: {outfile}") \ No newline at end of file diff --git a/inputs/Extrair_coords_time-lat-lon.gs b/inputs/Extract_coords_time-lat-lon.gs similarity index 100% rename from inputs/Extrair_coords_time-lat-lon.gs rename to inputs/Extract_coords_time-lat-lon.gs diff --git a/inputs/Extrair_coords_time-lat-lon_ciclone-automatico.gs b/inputs/Extract_coords_time-lat-lon_ciclone-automatico.gs similarity index 100% rename from inputs/Extrair_coords_time-lat-lon_ciclone-automatico.gs rename to inputs/Extract_coords_time-lat-lon_ciclone-automatico.gs diff --git a/inputs/namelist b/inputs/namelist index 3c04b75..9ddef92 100644 --- a/inputs/namelist +++ b/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 diff --git a/inputs/namelist_ERA5 b/inputs/namelist_ERA5 index d3624f8..9ddef92 100644 --- a/inputs/namelist_ERA5 +++ b/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 diff --git a/src/cli_interface.py b/src/cli_interface.py index 06e2db4..23daa4a 100644 --- a/src/cli_interface.py +++ b/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 \ No newline at end of file + return args diff --git a/src/data_handling.py b/src/data_handling.py index 0410c9d..4ed4b98 100644 --- a/src/data_handling.py +++ b/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 diff --git a/src/select_domain.py b/src/select_domain.py index 9c7aae5..2a40cad 100644 --- a/src/select_domain.py +++ b/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 \ No newline at end of file + return current_domain_limits diff --git a/src/utils.py b/src/utils.py index fcd90c9..92b26ba 100644 --- a/src/utils.py +++ b/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 \ No newline at end of file + return min_max_zeta, min_max_hgt, max_wind diff --git a/src/visualization.py b/src/visualization.py index 4af643e..116f99f 100644 --- a/src/visualization.py +++ b/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}") \ No newline at end of file + app_logger.error(f"❌ Failed to create Hovmöller diagram of mean zeta: {e}")