Scripts for WW3 data pre/post processing

[yunfangsun]

@aliabdolali shared a private repo (https://github.com/erdc/wave-tools) containing pre/post scripts for spectral wave models.

@saeed-moghimi-noaa

Attached is one example of directional spec interpolation from ascii to netcdf.

clear all
addpath /Users/rdchlaa9/Desktop/Denise_Grid/wave-tools/bin/matlab/
% read asii spec
[IN] = swden_ww3_read_ascii('bnd_ndbc.spec.spc');
%% user defined interpolation parameters
% input data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nfreq=32; % number of frequencies
nDir=36; % number of Directions
inc=1.1; % frequency increment
f0=0.050;
filenameIN='bnd_ndbc.spec_ww3.nc';% name of netcdf file (boundary)
filenameOUT='bnd_ndbc.spec_interp_ww3.nc';% name of netcdf file (boundary)
testcase='JXFL57'; % test vase
pointID='JXFL57'; % id (length should be 16)
coordinate='spherical'; % coordinate : Spherical, cartesian
visualize='true';
%----------------------------------------------------------%%----------------------------------------------------------%
deltatheta=360/nDir; % DeltaDir
%----------------------------------------------------------%
% convert ascii spec to nc spec
[filename] = write_directional_spectra_nc(filenameIN,testcase,...
IN.pointID,IN.lat,IN.lon,IN.dpt,IN.wnd,IN.wnddir,IN.cur,IN.curdir,...
IN.time,IN.frequency,IN.direction,IN.efth,coordinate);

% interpolate input spec to target spectral resolution
[IN,OUT]= spec_interp_nc(filenameIN,filenameOUT,nfreq,f0,inc,nDir,...
coordinate,visualize);

This script is using swden_ww3_read_ascii:

function [SWDEN] = swden_ww3_read_ascii(specfile)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function reads WW3 directional spectral density file %
% .spec asill format                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      Ali Abdolali March 2024 ali.abdolali@usace.army.mi   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  input data %--------------------------------------------%
% specfile: name of .spec file
%  output data %--------------------------------------------%
% buoy_name: list of buoy
% time (Matlab time)
% longitude of points (Degree)
% latitude of points (Degree)
% f: frequency (Hz)
% Dir: direction (degree)
% EFTH: Directional Spectral Density [freq,direction,stations,time]
% SPEC: Spectral Density [freq,stations,time]
% Hs: Significant Wave Heigth (m)
% Fp: Peak Freq (Hz)
% mwvdir: mean wave direction (degrees)
% cur: current speed (m/s)
% curdir: current direction (deg)
% wnd: wind speed (m/s)
% wnddir: wind direction (deg)
% dpt: depth (m)

%----------------------------------------------------------%
m=0;
%read and put all lines into cells
fid = fopen(specfile,'rt');
tline = fgetl(fid);
while ischar(tline)
m=m+1;
     file{m}=(tline);
    tline = fgetl(fid);
end
fclose(fid);
%----------------------------------------------------------%
freq=[];
Dir=[];
%'---------------------'    nfreq nDir   nstation '------------------------'
%'WAVEWATCH III SPECTRA'     50    36     1 'spectral resolution for points'
C = textscan(file{1},'%s %s %s %n %n %n %s %s %s %s');
nfreq=C{4};
nDir=C{5};
nStation=C{6};
  i=1;
    while length(freq)<nfreq
        i=i+1;
        A=cell2mat(textscan(file{i}, '%n', 'delimiter', '\n', 'whitespace', ' '));
    freq=[freq;A];
    end
%----------------------------------------------------------%    
  while length(Dir)<nDir
        i=i+1;
        A=cell2mat(textscan(file{i}, '%n', 'delimiter', '\n', 'whitespace', ' '));
    Dir=[Dir;A];
  end
%----------------------------------------------------------%
%go through time loop
  m=1;
   while i<length(file)
  i=i+1;
      time(m,1) = datenum(file{i},'yyyymmdd HHMMSS');
   i=i+1;
% buoyname  '  lat    long       dpt     wnd  wnddir cur  curdir   
%'46069     '  33.65-120.20     919.8   5.10  30.0   0.38 356.0
   % C=textscan(file{i}, '%s %s %n %n %n %n %n %n %n %*[^\n]');
   % buoy_name=C{1};
   % lat(1,m)=C{3};
   % lon(1,m)=C{4};
   % dpt(1,m)=C{5};
   % wnd(1,m)=C{6};
   % wnddir(1,m)=C{7};
   % cur(1,m)=C{8};
   % curdir(1,m)=C{9};
 C=textscan(file{9}(12+1:end), '%n %n %n %n %n %n %n %*[^\n]');
 buoy_name=sscanf(file{9}(2:11), '%s');
    lat(1,m)=C{1};
    lon(1,m)=C{2};
    dpt(1,m)=C{3};
    wnd(1,m)=C{4};
    wnddir(1,m)=C{5};
    cur(1,m)=C{6};
    curdir(1,m)=C{7};
%----------------------------------------------------------%    
  dens=[];
    while length(dens)<nDir*nfreq
        i=i+1;
        A=cell2mat(textscan(file{i}, '%n', 'delimiter', '\n', 'whitespace', ''));
    dens=[dens; A];
    end
  denstmp(1,:)=dens;
  EFTH(:,:,nStation,m)=fliplr(reshape(dens,[nDir,nfreq]));
   m=m+1;
   end
  
   
   
SWDEN.buoy_name=cellstr(buoy_name);
%SWDEN.pointID=cellstr(buoy_name);
SWDEN.pointID=[buoy_name,repmat(' ', [1, 16-strlength(buoy_name)])];
SWDEN.time=time;
SWDEN.lat=lat;
SWDEN.lon=lon;
SWDEN.frequency=freq;
SWDEN.direction=round(Dir*180/pi);%precision is not enough in spec
SWDEN.efth=EFTH;
SWDEN.cur=cur;
SWDEN.curdir=curdir;
SWDEN.wnd=wnd;
SWDEN.wnddir=wnddir;
SWDEN.dpt=dpt;
DeltaDir=abs(SWDEN.direction(2)-SWDEN.direction(1)); 
%DeltaDir
EFTH=SWDEN.efth;
EFTH=[EFTH; EFTH(1,:,:,:)];
[DD,FF,KK,TT]=size(SWDEN.efth);

EFTHCOS=cosd([SWDEN.direction;SWDEN.direction(1)]).*EFTH;
EFTHSIN=sind([SWDEN.direction;SWDEN.direction(1)]).*EFTH;

for k=1:KK
for ii=1:FF
SPEC(ii,k,:) = abs(trapz(0:pi*DeltaDir/180:2*pi,EFTH(:,ii,k,:)));
SPECA(ii,k,:) = (trapz(0:pi*DeltaDir/180:2*pi,EFTHCOS(:,ii,k,:)));
SPECB(ii,k,:) = (trapz(0:pi*DeltaDir/180:2*pi,EFTHSIN(:,ii,k,:)));
end
end
SWDEN.ef=SPEC;
for k=1:KK
 Hs(k,:)=4.004*sqrt(abs(trapz(2*pi*SWDEN.frequency,SPEC(:,k,:))))/sqrt(2*pi);
 AA(k,:)=((trapz(2*pi*SWDEN.frequency,SPECA(:,k,:))));
 BB(k,:)=((trapz(2*pi*SWDEN.frequency,SPECB(:,k,:))));
end
 SWDEN.hs=Hs;
 SWDEN.th1m=(180*atan2(BB,AA)/pi)+180;

for k=1:KK
    clear SPECmax
    clear SPECtmp
    SPECtmp(:,:)=SPEC(:,k,:);
    SPECmax=nanmax(SPECtmp);
for i=1:length(SWDEN.time)
    clear i1
    clear j1
    [i1,j1]=find(SPECtmp(:,i)==SPECmax(i));
    Fp(k,i)=nanmean(SWDEN.frequency(i1));
end
SWDEN.fp=Fp;    
end

write_directional_spectra_nc is as follow:

function [ncfile] = write_directional_spectra_nc(ncfile,testcase,...
    pointID,Lat,Lon,dpt,wndspd,wnddir,curspd,curdir,time,freq,dir,EFTH,coordinate)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function writes Directional Spectral density file in %
% WW3 netcdf format                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      Ali Abdolali April 2017 ali.abdolali@noaa.gov        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% testcase name: i.e. Inlet test'
% time
% ntime: length of time
% nfreq: number of frequency band
% nDir: number of directional band
% pointnumber : number of points- use 1 for this variable
% freq: frequency (Hz)
% dir: direction (rad)
% pointID: the boundary point name: i.e. 'b42001' 
% Lat: latitude (degrees)
% Lon: longitude (degrees)
% dpt: depth (m) [pointnumber x 1]
% wndspd: wind  speed (m/s) [pointnumber, ntime]
% wnddir: wind direction (degrees) [pointnumber, ntime]
% curspd: current speed (m/s) [pointnumber, ntime]
% curdir: current direction (degrees) [pointnumber, ntime]
% EFTH: directional spectral density [nDir x nfreq  x pointnumber x ntime]
% coordinate: options: 'spherical', 'cartesian'
% ncfile: name of output file 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%example: [ncfile] = write_directional_spectra_nc('B42001.nc',...
%'inlet','B42001',42,221.1,1521,30,270,0.5,36,time,freq,dir,EFTH)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nDir,nfreq,nStation,ntime] = size(EFTH);
frequency1=freq*1.047619000313776;
frequency1(1)=1;
frequency2=freq*0.952380928728317;
frequency2(end)=1;

nc = netcdf.create(ncfile, '64BIT_OFFSET');

% define global attributes
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'start_date',datestr(time(1)))
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'stop_date',datestr(time(end)))
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'source',testcase)
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'field_type','n/a')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'data_type','OCO spectra 2D')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'southernmost_latitude','n/a')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'northernmost_latitude','n/a')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'latitude_resolution','n/a')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'westernmost_longitude','n/a')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'easternmost_longitude','n/a')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'longitude_resolution','n/a')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'minimum_altitude','n/a')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'maximum_altitude','n/a')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'altitude_resolution','n/a')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'area','Global')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'format_version','1.1')
                  

% define dimensions

unlim_id = netcdf.defDim(nc, 'time', netcdf.getConstant('NC_UNLIMITED'));
STAT = netcdf.defDim(nc, 'station', nStation);
STRING = netcdf.defDim(nc, 'string16', 16);
FREQ = netcdf.defDim(nc, 'frequency', nfreq);
DIR = netcdf.defDim(nc, 'direction', nDir);

STATION=pointID;

time_varid = netcdf.defVar(nc, 'time', 'double', unlim_id);
netcdf.putAtt(nc, time_varid, 'long_name', 'julian day (UT)');
netcdf.putAtt(nc, time_varid, 'units', 'days since 1990-01-01 00:00:00');
netcdf.putAtt(nc, time_varid, 'field', 'time');
netcdf.putAtt(nc, time_varid, 'conventions', 'Relative julian days with decimal part (as parts of the day )');
netcdf.putAtt(nc, time_varid, 'axis', 'T');

station_varid = netcdf.defVar(nc, 'station', 'NC_INT', [STAT]);
netcdf.putAtt(nc, station_varid, 'long_name', 'station id');
netcdf.putAtt(nc, station_varid, 'axis', 'X');

station16_varid = netcdf.defVar(nc, 'string16', 'NC_INT', [STRING]);
netcdf.putAtt(nc, station_varid, 'long_name', 'station_name number of characters');
netcdf.putAtt(nc, station_varid, 'axis', 'W');

station_name_varid = netcdf.defVar(nc, 'station_name', 'NC_CHAR', [STAT, STRING]);
netcdf.putAtt(nc, station_varid, 'long_name', 'station name');
netcdf.putAtt(nc, station_varid, 'axis', 'XW');

tf=strcmp(coordinate,'spherical');                 

if tf==1
lon_varid=netcdf.defVar(nc, 'longitude' ,'NC_FLOAT',[STAT,unlim_id]);
netcdf.putAtt(nc, lon_varid, 'long_name', 'longitude');
netcdf.putAtt(nc, lon_varid, 'standard_name', 'longitude');
netcdf.putAtt(nc, lon_varid, 'globwave_name', 'longitude');
netcdf.putAtt(nc, lon_varid, 'content', 'TX');
netcdf.putAtt(nc, lon_varid, 'associates', 'time station');
netcdf.putAtt(nc, lon_varid, 'units', 'degree_east');

lat_varid=netcdf.defVar(nc, 'latitude' ,'NC_FLOAT',[STAT,unlim_id]);
netcdf.putAtt(nc, lat_varid, 'long_name', 'latitude');
netcdf.putAtt(nc, lat_varid, 'standard_name', 'latitude');
netcdf.putAtt(nc, lat_varid, 'globwave_name', 'latitude');
netcdf.putAtt(nc, lat_varid, 'content', 'TX');
netcdf.putAtt(nc, lat_varid, 'associates', 'time station');
netcdf.putAtt(nc, lat_varid, 'units', 'degree_north');
end

tf=strcmp(coordinate,'cartesian');                 
if tf==1
lon_varid=netcdf.defVar(nc, 'x' ,'NC_FLOAT',[STAT,unlim_id]);
netcdf.putAtt(nc, lon_varid, 'long_name', 'x');
netcdf.putAtt(nc, lon_varid, 'standard_name', 'x');
netcdf.putAtt(nc, lon_varid, 'globwave_name', 'x');
netcdf.putAtt(nc, lon_varid, 'content', 'TX');
netcdf.putAtt(nc, lon_varid, 'associates', 'time station');
netcdf.putAtt(nc, lon_varid, 'units', 'm');

lat_varid=netcdf.defVar(nc, 'y' ,'NC_FLOAT',[STAT,unlim_id]);
netcdf.putAtt(nc, lat_varid, 'long_name', 'y');
netcdf.putAtt(nc, lat_varid, 'standard_name', 'y');
netcdf.putAtt(nc, lat_varid, 'globwave_name', 'y');
netcdf.putAtt(nc, lat_varid, 'content', 'TX');
netcdf.putAtt(nc, lat_varid, 'associates', 'time station');
netcdf.putAtt(nc, lat_varid, 'units', 'm');
end


f_varid=netcdf.defVar(nc, 'frequency' ,'NC_FLOAT',[FREQ]);
netcdf.putAtt(nc, f_varid, 'long_name', 'frequency of center band');
netcdf.putAtt(nc, f_varid, 'standard_name', 'sea_surface_wave_frequency');
netcdf.putAtt(nc, f_varid, 'globwave_name', 'frequency');
netcdf.putAtt(nc, f_varid, 'units', 's-1');
netcdf.putAtt(nc, f_varid, 'axis', 'Y');

f1_varid=netcdf.defVar(nc, 'frequency1' ,'NC_FLOAT',[FREQ]);
netcdf.putAtt(nc, f1_varid, 'long_name', 'frequency of lower band');
netcdf.putAtt(nc, f1_varid, 'standard_name', 'frequency_of_lower_band');
netcdf.putAtt(nc, f1_varid, 'globwave_name', 'frequency_lower_band');
netcdf.putAtt(nc, f1_varid, 'units', 's-1');
netcdf.putAtt(nc, f1_varid, 'axis', 'Y');
netcdf.putAtt(nc, f1_varid, 'associates', 'frequency');

f2_varid=netcdf.defVar(nc, 'frequency2' ,'NC_FLOAT',[FREQ]);
netcdf.putAtt(nc, f2_varid, 'long_name', 'frequency of upper band');
netcdf.putAtt(nc, f2_varid, 'standard_name', 'frequency_of_upper_band');
netcdf.putAtt(nc, f2_varid, 'globwave_name', 'frequency_upper_band');
netcdf.putAtt(nc, f2_varid, 'units', 's-1');
netcdf.putAtt(nc, f2_varid, 'axis', 'Y');
netcdf.putAtt(nc, f2_varid, 'associates', 'frequency');


dir_varid=netcdf.defVar(nc, 'direction' ,'NC_FLOAT',[DIR]);
netcdf.putAtt(nc, dir_varid, 'long_name', 'sea surface wave to direction');
netcdf.putAtt(nc, dir_varid, 'standard_name', 'sea surface wave to direction');
netcdf.putAtt(nc, dir_varid, 'globwave_name', 'direction');
netcdf.putAtt(nc, dir_varid, 'units', 'degree');
netcdf.putAtt(nc, time_varid, 'axis', 'X');

EFTH_varid=netcdf.defVar(nc, 'efth' ,'NC_FLOAT',[DIR,FREQ,STAT,unlim_id]);
netcdf.putAtt(nc, EFTH_varid, 'units', 'm2 s rad-1');
netcdf.putAtt(nc, EFTH_varid, 'field', 'efth, scalar, series');
netcdf.putAtt(nc, EFTH_varid, 'long_name', 'sea surface wave directional variance spectral density');
netcdf.putAtt(nc, EFTH_varid, 'standard_name', 'sea_surface_wave_directional_variance_spectral_density');
netcdf.putAtt(nc, EFTH_varid, 'globwave_name', 'directional_variance_spectral_density');
netcdf.putAtt(nc, EFTH_varid, 'associates', 'time station frequency direction');
netcdf.putAtt(nc, EFTH_varid, 'content', 'TXYZ');


d_varid=netcdf.defVar(nc, 'dpt' ,'NC_FLOAT',[STAT,unlim_id]);
netcdf.putAtt(nc, d_varid, 'long_name', 'depth');
netcdf.putAtt(nc, d_varid, 'standard_name', 'depth');
netcdf.putAtt(nc, d_varid, 'globwave_name', 'depth');
netcdf.putAtt(nc, d_varid, 'content', 'TX');
netcdf.putAtt(nc, d_varid, 'associates', 'time station');
netcdf.putAtt(nc, d_varid, 'units', 'm');


wnd_varid=netcdf.defVar(nc, 'wnd' ,'NC_FLOAT',[STAT,unlim_id]);
netcdf.putAtt(nc, wnd_varid, 'units', 'm s-1');
netcdf.putAtt(nc, wnd_varid, 'long_name', 'wind speed at 10m');
netcdf.putAtt(nc, wnd_varid, 'standard_name', 'wind speed');
netcdf.putAtt(nc, wnd_varid, 'globwave_name', 'wind speed');
netcdf.putAtt(nc, wnd_varid, 'content', 'TX');
netcdf.putAtt(nc, wnd_varid, 'associates', 'time station');


wnddir_varid=netcdf.defVar(nc, 'wnddir' ,'NC_FLOAT',[STAT,unlim_id]);
netcdf.putAtt(nc, wnddir_varid, 'units', 'degree');
netcdf.putAtt(nc, wnddir_varid, 'long_name', 'wind direction');
netcdf.putAtt(nc, wnddir_varid, 'standard_name', 'wind_from_direction');
netcdf.putAtt(nc, wnddir_varid, 'globwave_name', 'wind_from_direction');
netcdf.putAtt(nc, wnddir_varid, 'content', 'TX');
netcdf.putAtt(nc, wnddir_varid, 'associates', 'time station');


cur_varid=netcdf.defVar(nc, 'cur' ,'NC_FLOAT',[STAT,unlim_id]);
netcdf.putAtt(nc, cur_varid, 'units', 'm s-1');
netcdf.putAtt(nc, cur_varid, 'long_name', 'sea water speed');
netcdf.putAtt(nc, cur_varid, 'standard_name', 'sea_water_speed');
netcdf.putAtt(nc, cur_varid, 'globwave_name', 'sea_water_speed');
netcdf.putAtt(nc, cur_varid, 'content', 'TX');
netcdf.putAtt(nc, cur_varid, 'associates', 'time station');


curdir_varid=netcdf.defVar(nc, 'curdir' ,'NC_FLOAT',[STAT,unlim_id]);
netcdf.putAtt(nc, curdir_varid, 'units', 'degree');
netcdf.putAtt(nc, curdir_varid, 'long_name', 'direction from of sea water velocity');
netcdf.putAtt(nc, curdir_varid, 'standard_name', 'direction_of_sea_water_velocity');
netcdf.putAtt(nc, curdir_varid, 'globwave_name', 'direction_of_sea_water_velocity');
netcdf.putAtt(nc, curdir_varid, 'content', 'TX');
netcdf.putAtt(nc, curdir_varid, 'associates', 'time station');
                     


netcdf.endDef(nc);

netcdf.putVar(nc, time_varid,0,length(time), time-datenum('01011990 000000','ddmmyyyy HHMMSS'));
netcdf.putVar(nc, station_varid, 84);
netcdf.putVar(nc, station16_varid, nan*(ones(length(pointID),1)));
netcdf.putVar(nc, station_name_varid,pointID);
netcdf.putVar(nc, lon_varid, Lon);
netcdf.putVar(nc, lat_varid, Lat);
netcdf.putVar(nc, f_varid, freq);
netcdf.putVar(nc, f1_varid, frequency1);
netcdf.putVar(nc, f2_varid, frequency2);
netcdf.putVar(nc, dir_varid, dir);
netcdf.putVar(nc, EFTH_varid, EFTH);
netcdf.putVar(nc, d_varid, dpt);
netcdf.putVar(nc, wnd_varid, wndspd);
netcdf.putVar(nc, wnddir_varid, wnddir);
netcdf.putVar(nc, cur_varid, curspd);
netcdf.putVar(nc, curdir_varid, curdir);




netcdf.close(nc);


%%

fileattrib(ncfile,'+w');
ncwriteatt(ncfile,'efth','scale_factor', 1);
ncwriteatt(ncfile,'efth','add_offset', 0);
ncwriteatt(ncfile,'efth','valid_min', 0);
ncwriteatt(ncfile,'efth','valid_max', 1e+20);
ncwriteatt(ncfile,'efth','_FillValue', int32(9.97e+36));

ncwriteatt(ncfile,'dpt','scale_factor', 1);
ncwriteatt(ncfile,'dpt','add_offset', 0);
ncwriteatt(ncfile,'dpt','valid_min', -100);
ncwriteatt(ncfile,'dpt','valid_max', 1e+04);
ncwriteatt(ncfile,'dpt','_FillValue', int32(9.97e+36));

tf=strcmp(coordinate,'spherical');              
if tf==1
ncwriteatt(ncfile,'latitude','scale_factor', 1);
ncwriteatt(ncfile,'latitude','add_offset', 0);
ncwriteatt(ncfile,'latitude','valid_min', -90);
ncwriteatt(ncfile,'latitude','valid_max', 180);
ncwriteatt(ncfile,'latitude','_FillValue', int32(9.97e+36));

ncwriteatt(ncfile,'longitude','scale_factor', 1);
ncwriteatt(ncfile,'longitude','add_offset', 0);
ncwriteatt(ncfile,'longitude','valid_min', -180);
ncwriteatt(ncfile,'longitude','valid_max', 360);
ncwriteatt(ncfile,'longitude','_FillValue', int32(9.97e+36));
end

tf=strcmp(coordinate,'cartesian');
if tf==1
ncwriteatt(ncfile,'y','scale_factor', 1);
ncwriteatt(ncfile,'y','add_offset', 0);
ncwriteatt(ncfile,'y','valid_min', -90);
ncwriteatt(ncfile,'y','valid_max', 180);
ncwriteatt(ncfile,'y','_FillValue', int32(9.97e+36));

ncwriteatt(ncfile,'x','scale_factor', 1);
ncwriteatt(ncfile,'x','add_offset', 0);
ncwriteatt(ncfile,'x','valid_min', -180);
ncwriteatt(ncfile,'x','valid_max', 360);
ncwriteatt(ncfile,'x','_FillValue', int32(9.97e+36));
end


ncwriteatt(ncfile,'direction','scale_factor', 1);
ncwriteatt(ncfile,'direction','add_offset', 0);
ncwriteatt(ncfile,'direction','valid_min', 0);
ncwriteatt(ncfile,'direction','valid_max', 360);
ncwriteatt(ncfile,'direction','_FillValue', int32(9.97e+36));

                          
ncwriteatt(ncfile,'frequency','scale_factor', 1);
ncwriteatt(ncfile,'frequency','add_offset', 0);
ncwriteatt(ncfile,'frequency','valid_min', 0);
ncwriteatt(ncfile,'frequency','valid_max', 10);
ncwriteatt(ncfile,'frequency','_FillValue', int32(9.97e+36));

    
ncwriteatt(ncfile,'frequency1','scale_factor', 1);
ncwriteatt(ncfile,'frequency1','add_offset', 0);
ncwriteatt(ncfile,'frequency1','valid_min', 0);
ncwriteatt(ncfile,'frequency1','valid_max', 10);
ncwriteatt(ncfile,'frequency1','_FillValue', int32(9.97e+36));


ncwriteatt(ncfile,'frequency2','scale_factor', 1);
ncwriteatt(ncfile,'frequency2','add_offset', 0);
ncwriteatt(ncfile,'frequency2','valid_min', 0);
ncwriteatt(ncfile,'frequency2','valid_max', 10);
ncwriteatt(ncfile,'frequency2','_FillValue', int32(9.97e+36));


ncwriteatt(ncfile,'wnd','scale_factor', 1);
ncwriteatt(ncfile,'wnd','add_offset', 0);
ncwriteatt(ncfile,'wnd','valid_min', 0);
ncwriteatt(ncfile,'wnd','valid_max', 100);
ncwriteatt(ncfile,'wnd','_FillValue', int32(9.97e+36));

ncwriteatt(ncfile,'wnddir','scale_factor', 1);
ncwriteatt(ncfile,'wnddir','add_offset', 0);
ncwriteatt(ncfile,'wnddir','valid_min', 0);
ncwriteatt(ncfile,'wnddir','valid_max', 360);
ncwriteatt(ncfile,'wnddir','_FillValue', int32(9.97e+36));


ncwriteatt(ncfile,'cur','scale_factor', 1);
ncwriteatt(ncfile,'cur','add_offset', 0);
ncwriteatt(ncfile,'cur','valid_min', 0);
ncwriteatt(ncfile,'cur','valid_max', 100);
ncwriteatt(ncfile,'cur','_FillValue', int32(9.97e+36));

ncwriteatt(ncfile,'curdir','scale_factor', 1);
ncwriteatt(ncfile,'curdir','add_offset', 0);
ncwriteatt(ncfile,'curdir','valid_min', 0);
ncwriteatt(ncfile,'curdir','valid_max', 360);
ncwriteatt(ncfile,'curdir','_FillValue', int32(9.97e+36));
ncwriteatt(ncfile,'station','_FillValue', int32(-2147483647));
ncwriteatt(ncfile,'string16','_FillValue', int32(-2147483647));
end

the interpolation
function[IN,OUT]=spec_interp_nc(filenameIN,filenameOUT,nfreq,f0,inc,nDir,coordinate,visualize)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program reads all netcdf file contents %
% ali.abdolali@erdc.dren.mil March 06, 2024 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% INPUT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ncf: name of input netcdf file
% nfreq=35; % number of frequencies
% nDir=36; % number of Directions
% inc=1.1; % frequency increment
% f0=0.038; % first freq
% visualize='true';
% filenameIN='JXFL57_swan.nc';% name of netcdf file (boundary)
% coordinate='spherical'; % coordinate : Spherical, cartesian
%%%%%%%%%%%%%%%%%%% OUTPUT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% filenameOUT='JXFL57_ww3.nc';% name of netcdf file (boundary)
%%%%%%%%%%%%%%%%%%%%%%%% Dependency %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert_time
% read_nc
% dictionary_wave
%%%%%%%%%%%%%%%%%%% example %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[IN,OUT]= spec_interp_nc('in.nc','out.nc',35,0.038,1.1,36,...
% 'spherical','true')
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------%
testcase=['interpolated from ',filenameIN]; % test vase
deltatheta=360/nDir; % DeltaDir
% Dir0=5; % first direction (deg)
% theta0=0; % first Dir
%-------------------------------------------------------------------------%
% read input
IN=read_nc(filenameIN);
% sort direction and rearrange efth
[IN.direction,id]=sort(IN.direction);
IN.efth=IN.efth(id,:,:);
%-------------------------------------------------------------------------%
% prepare freq and direction
f=f0;
for i=1:nfreq-1
f(end+1)=f(end)inc;
end
dir0(1,:)=90:-deltatheta:0;
dir0(end+1:nDir)=360-deltatheta+dir0(end):-deltatheta:90+deltatheta-dir0(end);
dir(1,:)=pidir0/180;
[dirtmpO,ftmpO]=meshgrid(dir0,f);
%-------------------------------------------------------------------------%
% create output structure
OUT=IN;
OUT=rmfield(OUT,'efth');
OUT.pointID=[[OUT.station_name{:}],repmat(' ', [1, 16-strlength([OUT.station_name{:}])])];
%-------------------------------------------------------------------------%
% make a tmp IN and add two columns on both sides of directions
INtmp=IN;
INtmp=rmfield(INtmp,'direction');
INtmp=rmfield(INtmp,'efth');
imax=find(IN.direction==max(IN.direction(:)));
imin=find(IN.direction==min(IN.direction(:)));
INtmp.direction=[IN.direction(imax)-360; IN.direction;IN.direction(imin)+360];
INtmp.efth=[IN.efth(imax,:,:); IN.efth; IN.efth(imin,:,:)];
%-------------------------------------------------------------------------%
% interpolation
[dirtmp,ftmp]=meshgrid(INtmp.direction,INtmp.frequency);
for i=1:length(IN.time)
DENStmp(:,:)=INtmp.efth(:,:,i);
DENStmp2=DENStmp';
tmp=reshape(interp2(dirtmp,ftmp,DENStmp',dirtmpO(:),ftmpO(:)),size(dirtmpO));
OUT.efth(:,:,1,i)=tmp';
end
%-------------------------------------------------------------------------%
display (['Generating ', filenameOUT,' ...'])
%dump into netcdf
[filename] = write_directional_spectra_nc(filenameOUT,testcase,...
OUT.pointID,OUT.lat,OUT.lon,OUT.dpt,OUT.wnd,OUT.wnddir,...
OUT.cur,OUT.curdir,OUT.time,f,dir0,OUT.efth,coordinate);
%%
%%
tf=strcmp(visualize,'true');
if tf==1
display (['Plot Generation ', [filename,'.gif'],' ...'])
width=500; % Width of figure for movie [pixels]
height=800; % Height of figure of movie [pixels]
left=200; % Left margin between figure and screen edge [pixels]
bottom=200; % Bottom margin between figure and screen edge [pixels]

h=figure('Color','w');
set(gcf,'Position', [left bottom width height])
for k=1:length(OUT.time)
k
clf

sss1=subplot(2,1,1);
clear SPEC1
clear SPECC
SPEC1(:,:)=IN.efth(:,:,k);
SPECC=SPEC1;
[~,cc] = polarPcolor(IN.frequency',[IN.direction; IN.direction(1)]',...
[SPECC; SPECC(1,:)],'Nspokes',36,'ncolor',10,'labelR','f (Hz)','Rscale','log');
ylabel(cc,['INPUT - ',datestr(IN.time(k))],'FontSize',14);

sss2=subplot(2,1,2);
clear SPECC
SPECC=OUT.efth(:,:,1,k);
[~,cc] = polarPcolor(f,[dir0 dir0(1)],...
[SPECC; SPECC(1,:)],'Nspokes',36,'ncolor',10,'labelR','f (Hz)','Rscale','log');
ylabel(cc,['OUTPUT - ',datestr(OUT.time(k))],'FontSize',14);

set(sss1, 'Position', [.08 .55 .88 .4]);
set(sss2, 'Position', [.08 .08 .88 .4]);
frame = getframe(h);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);

if k == 1
    imwrite(imind,cm,[filename,'.gif'],'gif', 'Loopcount',inf);
else
    imwrite(imind,cm,[filename,'.gif'],'gif','WriteMode','append');
end   

end
end
%----------------------------------------------------------%
display (['Finished'])

[saeed-moghimi-noaa]

Hi @yunfangsun @aliabdolali

Please continue the discussion here.

Thanks,
-Saeeed

[yunfangsun]

@aliabdolali, I can run your example case and reproduce the results, however, when I use the GFS wave spec which I am using in my configuration, it produces the following error for the ascii read ([IN] = swden_ww3_read_ascii('../obc/gfswave.22108.spec');):

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in swden_ww3_read_ascii (line 87)
    curdir(1,m)=C{7};

Error in convert_WW3ascii_interp_2_WW3 (line 4)
[IN] = swden_ww3_read_ascii('../obc/gfswave.22108.spec');

could you please help me fix the error for those files?

[aliabdolali]

please pull wave-tools and try again.

[yunfangsun]

Hi @aliabdolali ,

Thank you, that error disappeared, however, the results seem there is a mismatch in the reading, please take a look at the figure,

[aliabdolali]

You need to be careful about spectral resolution, WW3 freq increment is logarithmic (the larger the frequeny, the wider the spectral band), and in GFS, it starts with 0.035 with 1.07.
If you want to interpolate it to 32 direction with coarser increment (1.1), then you will loose some thing due to coarsening.
BTW, do not use what I used for spectral resoltion, you need to define your target spectral resolution:
nfreq=32?; % number of frequencies
nDir=36?; % number of Directions
inc=1.1; % frequency increment
f0=0.0350?;
check the bold ones with your target.
Let me know if you need further explanation.

[aliabdolali]

when you read netcdf files, it automatically calculates hs, tp, ..., so you can compare them to make sure the mismatch is acceptable or not.
The mismatch is not sth that we can fix, when you coarsen your spectral resolution, you will lose accuracy.

[yunfangsun]

#8

[janahaddad]

will continue updating here as needed as @yunfangsun builds scripts/notebooks for WW3 model pre- and post-processing.

[janahaddad]

Update from Yunfang: updating scripts to prep namelist generation