Open yunfangsun opened 3 months ago
Hi @yunfangsun @aliabdolali
Please continue the discussion here.
Thanks, -Saeeed
@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?
please pull wave-tools and try again.
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,
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.
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.
will continue updating here as needed as @yunfangsun builds scripts/notebooks for WW3 model pre- and post-processing.
Update from Yunfang: updating scripts to prep namelist generation
@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:
write_directional_spectra_nc is as follow:
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);
end end %----------------------------------------------------------%
display (['Finished'])