UFS-CAT Training: Check RT SCHISM output

[janahaddad]

Problem

This issue is focused on prepping for the UFS-CAT training session. The currently passing RTs (DATM+SCHISM, DATM+WW3) have only been tested by Ufuk, and need to be re-run by another team member who can plot and verify the model outputs.

Solution

re-run coastal_ike_shinnecock_atm2sch on Hercules: RT directory: /work2/noaa/nems/tufuk/RT Baseline Directory: /work2/noaa/nems/tufuk/RT/NEMSfv3gfs/develop-20240126

plot SCHISM model output and verify results

[janahaddad]

related to oceanmodeling/schism#5

[yunfangsun]

To plot the water elevation field cartoon, the function is as follows:

def plot_ocn_from_file(file_path, vara, vmin, vmax, output_dir='output_images'):

Load the external mesh for triangulation

hgrid = pyschism.mesh.Hgrid.open('rt_tests_figures/rt_990727/coastal_ike_shinnecock_atm2sch2ww3_intel/hgrid.gr3', crs=4326)
triangulation = hgrid.triangulation

# Open the NetCDF file
ds = nc.Dataset(file_path)

# Get time variable and convert to string format for filenames
time_var = ds.variables['time']
time_units = time_var.units

# Extract long name and units for the variable

long_name = ds.variables[vara].getncattr('long_name') if 'long_name' in ds.variables[vara].ncattrs() else vara

units = ds.variables[vara].getncattr('units') if 'units' in ds.variables[vara].ncattrs() else 'no units'

# Loop through each time step
for t_idx in range(len(time_var)):
    el = ds.variables[vara][t_idx, :].flatten()

    # Convert time value to a string
    time_val = nc.num2date(time_var[t_idx], time_units)
    time_str = time_val.strftime('%Y%m%d_%H%M%S')

    # Ensure the output directory exists
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    # Plotting
    plt.figure(figsize=(10, 6))
    plt.tripcolor(triangulation, el, shading='flat', cmap=cmap, vmin=vmin, vmax=vmax)
    plt.colorbar(label=f'Elevation (m)', extend='both', aspect=20, shrink=0.8, orientation='vertical')

    title_text = 'Elevation' + f' - {time_str}'
    plt.title(title_text)
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')

    # Save the figure
    file_name = f'{vara}_{time_str}.png'
    output_path = os.path.join(output_dir, file_name)
    plt.savefig(output_path)
    plt.close()

# Close the NetCDF dataset
ds.close()

To use it

cmap = LinearSegmentedColormap.from_list('red_blue_colormap', ['blue', 'white', 'red'])

directory_path = 'rt_tests_figures/rt_990727/coastal_ike_shinnecock_atm2sch2ww3_intel/outputs/'

file_paths = glob.glob(os.path.join(directorypath, 'schout?.nc'))

Loop through each file and plot

for file_path in file_paths: plot_ocn_from_file(file_path,'elev',-0.6,0.6)

!convert outputimages/elev*.png output_images/elev.gif

[yunfangsun]

To plot the surface forcing cartoon:

Open the NetCDF file

dataset = nc.Dataset('rt_tests_figures/rt_562057/coastal_ike_shinnecock_atm2sch_intel/INPUT/wind_atm_fin_ch_time_vec_STR_fixed.nc')

Extract variables

latitude = dataset.variables['latitude'][:] longitude = dataset.variables['longitude'][:] time = dataset.variables['time'] time_units = time.units time_values = time[:] # Get all time values

Convert time values to datetime objects

dates = num2date(time_values, units=time_units)

Set up the figure and axis with cartopy

fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={'projection': ccrs.PlateCarree()})

Define the color normalization with fixed range 0 to 20 m/s (not 0 to 500 as in the original script)

norm = colors.Normalize(vmin=0, vmax=30)

Create a ScalarMappable for the colorbar with the fixed normalization

sm = plt.cm.ScalarMappable(cmap='coolwarm', norm=norm) sm.set_array([])

Add the colorbar to the figure

cbar = fig.colorbar(sm, ax=ax, orientation='vertical', extend='both', fraction=0.046, pad=0.04) cbar.set_label('Wind Speed (m/s)')

Function to draw each frame of the animation

def update(frame): ax.clear() # Clear previous frame ax.set_extent([-72.9, -72.04, 40.4, 40.99], crs=ccrs.PlateCarree())

# Re-add static features for each frame
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.RIVERS)

# Extract wind components and calculate speed
uwnd = dataset.variables['uwnd'][frame, :, :]
vwnd = dataset.variables['vwnd'][frame, :, :]
wind_speed = np.sqrt(uwnd**2 + vwnd**2)

# Plot wind speed as a contour map
speed_contour = ax.contourf(longitude, latitude, wind_speed, levels=np.linspace(0, 30, 31), cmap='coolwarm', transform=ccrs.PlateCarree())
# This line adds a colorbar for each frame, which can be redundant. Consider adding it outside the update function if necessary.
geom_poly_2 = geom_obj_2.get_multipolygon()

gpd.GeoSeries(geom_poly_2).plot(color='c', ax=ax)

plot_mesh_edge(mesh_obj_2, ax=ax, color='k', linewidth=0.1)
# Format the date to a readable string format, e.g., "2023-01-01 12:00"
date_str = dates[frame].strftime('%Y-%m-%d %H:%M')
plt.title(f'Wind Speed at {date_str}')

Create animation using the 'update' function

ani = animation.FuncAnimation(fig, update, frames=len(time_values), interval=200)

Save the animation

ani.save('rt_tests_figures/shinecock_atm/wind_speed_contour_map_corrected_ike.gif', writer='pillow', fps=5)

[janahaddad]

closing & tracking in https://github.com/oceanmodeling/ufs-coastal/issues/43 instead