UFS-CAT Training: Check RT WW3 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_atm2ww3 on Hercules: RT directory: /work2/noaa/nems/tufuk/RT Baseline Directory: /work2/noaa/nems/tufuk/RT/NEMSfv3gfs/develop-20240126

plot WW3 model output and verify results

[yunfangsun]

To make the wave filed from the hourly output, the function is as follows `def plot_ww_from_file(file_path, vara,vmin,vmax,output_dir='output_images'):

Open the netCDF file

ds = nc.Dataset(file_path)

# Extract variables
lon = ds.variables['lon'][:].data.flatten()
lat = ds.variables['lat'][:].data.flatten()
hs = ds.variables[vara][:].data.flatten()

# nv needs to be transposed if it's not in the (N, 3) shape
nv = ds.variables['nconn'][:] - 1  # Adjust for 0-based indexing if necessary
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'

ds.close()

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

# Create a Triangulation object
triangulation = tri.Triangulation(lon, lat, triangles=nv)
title_text = long_name + ' - ' + os.path.basename(file_path).replace('.nc', '')

# Plotting
plt.figure(figsize=(10, 6))
plt.tripcolor(triangulation, hs, shading='flat', cmap=cmap, vmin=vmin, vmax=vmax)
plt.colorbar(label=f'{long_name} ({units})', extend='both', aspect=20, shrink=0.8, orientation='vertical')

plt.title(title_text)
plt.xlabel('Longitude')
plt.ylabel('Latitude')

# Save the figure
prefix = vara  + '_'# The characters you want to add at the beginning
file_name = prefix + os.path.basename(file_path).replace('.nc', '.png')

output_path = os.path.join(output_dir, file_name)
plt.savefig(output_path)
plt.close()

`

To use the function:

import glob from matplotlib.colors import LinearSegmentedColormap

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

directory_path = 'rt_tests_figures/rt_990727/coastal_ike_shinnecock_atm2sch2ww3_intel/test'

file_paths = glob.glob(os.path.join(directory_path, '2008*.nc'))

Loop through each file and plot

for file_path in file_paths: plot_ww_from_file(file_path,'HS',0.1,0.25) plot_ww_from_file(file_path,'SXX',0,10) plot_ww_from_file(file_path,'SXY',0,10) plot_ww_from_file(file_path,'SYY',0,20)

!convert outputimages/HS.png output_images/hs.gif !convert outputimages/SXX.png output_images/sxx.gif !convert outputimages/SYY.png output_images/syy.gif !convert outputimages/SXY.png output_images/sxy.gif

[yunfangsun]

To plot the wave model-simulation comparisons:

def plot_wave_TimeSeries_comparison(data1=None, data2=None, data3=None, data_observed=None, **kwargs): rounding_base = kwargs.get('rounding_base', None) station_label = kwargs.get('station_label', None) DPI = kwargs.get('DPI', None) save_file_name = kwargs.get('save_file_name', None) date_rng = kwargs.get('date_rng', None) # in datetime format (datetime(yyyy, mm, dd))

# Figure dpi resolution:
if DPI is None:       DPI = 120

Nplots = 3 #data.shape[1]  # number of subplots (=number of selcted stations)

#
row, col = 1, Nplots # number of rows and columns
figW, figH = 16.5, min(8.0, row*2.00)  # Figure width and height
figWspace, figHsapce = 0.25, 0.15

fig, axs = plt.subplots(row , col, figsize=(figW, figH), facecolor='w', edgecolor='k')
fig.subplots_adjust(hspace = figHsapce, wspace=figWspace)

axs = axs.ravel()

for i in range(Nplots):

    if i == 0:
        if data1 is not None:
            axs[i].plot(data1.iloc[:], linewidth=0.5, label='model', color = 'red')  # predicted
    elif i == 1:
        if data2 is not None:
            axs[i].plot(data2.iloc[:], linewidth=0.5, label='model', color = 'red')             
    else:
        if data3 is not None:
            axs[i].plot(data3.iloc[:], linewidth=0.5, label='model', color = 'red')            

    if data_observed is not None:
        if i == 0:
            axs[i].plot(data_observed["wave_height"][date_rng[0]:date_rng[len(date_rng)-1]], linewidth=0.5, label='obs', color = 'black' )

        elif i == 2:
            arr = data_observed["mean_wave_direction"][date_rng[0]:date_rng[len(date_rng)-1]]
            tmp = arr.values
            tmp = np.where(tmp > 360.0, np.nan, tmp)
            arr.values[:] = tmp
            axs[i].plot(arr, linewidth=0.5, label='obs', color = 'black' )
        else:
            axs[i].plot(data_observed["dominant_period"][date_rng[0]:date_rng[len(date_rng)-1]], linewidth=0.5, label='obs', color = 'black' )

    axs[i].grid(axis = 'both', color = 'gray', linestyle = '-', linewidth = 0.75, alpha=0.15)
    axs[i].tick_params(axis="both",direction="in")  #, pad=0
    #plt.setp(axs[i].xaxis.get_majorticklabels(), rotation=90)   

    # beautify the x-labels
    plt.gcf().autofmt_xdate()
    # Define the date format
    date_form = DateFormatter("%b-%d")
    axs[i].xaxis.set_major_formatter(date_form)        
    axs[Nplots-2].tick_params(axis='x', rotation=90)
    axs[Nplots-1].tick_params(axis='x', rotation=90) 
    axs[Nplots-3].tick_params(axis='x', rotation=90)
    plt.setp(axs[Nplots-2].xaxis.get_majorticklabels(), ha='left')
    plt.setp(axs[Nplots-1].xaxis.get_majorticklabels(), ha='left')
    plt.setp(axs[Nplots-3].xaxis.get_majorticklabels(), ha='left')

    if date_rng is not None:
        plt.setp(axs[i], xticks=date_rng)
        axs[i].set_xlim([date_rng[0], date_rng[-1]])

    if i == 0:
        axs[i].set_ylim([0 , 12])
        axs[i].text(0.5, 0.975, 'SWH',
             horizontalalignment='center', verticalalignment='top',
             transform=axs[i].transAxes, size=7,weight='bold') 

    elif i == 2:
        axs[i].set_ylim([0, 360])
        #axs[i].legend()
        axs[i].text(0.5, 0.975, 'mean wave direction',
             horizontalalignment='center', verticalalignment='top',
             transform=axs[i].transAxes, size=7,weight='bold')              

    else:
        axs[i].set_ylim([0, 25])
        axs[i].text(0.5, 0.975, 'dominant_period',
             horizontalalignment='center', verticalalignment='top',
             transform=axs[i].transAxes, size=7,weight='bold')            

axs[Nplots-1].legend(loc="lower right", ncol=3)    
fig.add_subplot(111, frame_on=False)
plt.tick_params(labelcolor="none", bottom=False, left=False)

plt.title(station_label, size=11,weight='bold')

fig.savefig('wave.png', dpi=150, bbox_inches='tight', format='png')

[janahaddad]

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