make animation for fort.63.nc

[Jiangchao3]

Hi @krober10nd and @WPringle,

**Recently, I am trying to do some visualization for my simulation result.

I just completed an initial version for make animation for fort.63 file in nc format.

Here is the video and my code**

video.mp4

clearvars; clc;
addpath(genpath('/OceanMesh2D/'))
addpath(genpath('/OceanMesh2D/utilities/'))
Animation63nc('nc63','fort.63.nc', ...
'BestTrack','bio062007.txt', ...
'Base_date','2007-11-10 06:00:00', ...
'Title','Water Level Elevation During SIDR Hit Bangladesh', ...
'FigurePosition',[0.1,0.25,0.7,0.75], ...
'IterStart',1, ...
'IterStop',20, ...
'VideoWriteDir','OceanMesh2D/make_animation/')

function Animation63nc(varargin)
% Anim63nc generate an animatiom of an ADCIRC 63 netCDF file pair
% Some idears come from https://github.com/BrianOBlanton/adcirc_util
%
% P/V pairs:
% nc63 - Elevation Time Series at All Nodes in the Model Grid
% BestTrack -
% VideoWriteDir - Path to save the Video
% VideoProfile - mp4 use 'MPEG-4', but only in Windows and macOS
% AxisLims - plot axis limits (def=grid lims)
% Title - title string as a cell array (def={''})
% IterStart - starting iteration number (def=1)
% IterStride - iteration stride (skip) (def=1)
% IterStop - stopping iteration number (def=-1, represent end)
% ColorMin - min (def=min(zeta))
% ColorMax - max (def=max(zeta))
% ColorMap - colormap to use (def=jet(32))
% ImageReso - (def='-r200';)
% Base_date - datenum of starttime
% FigurePosition - we use normalized position
% Projection - see projections in m_map
% FontSize - default value is set as 16
%
% Author: Jiangchao Qiu, (MIT/ESSG; email:qiujch24@mit.edu)
% Created: May 14 2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p.nc63='';
p.BestTrack='';
p.VideoWriteDir='./';
p.VideoProfile='MPEG-4';
p.AxisLims=[];
p.Title={};
p.IterStart=1;
p.IterStride=1;
p.IterStop=-1;
p.ColorMin=[];
p.ColorMax=[];
p.ColorMap=jet(32);
p.ImageReso='-r200';
p.Base_date=[];
p.FigurePosition=[];
p.Projection='Miller';
p.FontSize=18;

p = parse_pv_pairs(p,varargin);

%% read into data from nc format file
x = ncread(p.nc63,'x');
y = ncread(p.nc63,'y');
element = ncread(p.nc63,'element')';
zeta = ncread(p.nc63,'zeta');
time = ncread(p.nc63,'time');
% find the max and min zeta value as the up and down limit of colormap
min_zeta = min(min(zeta)); p.ColorMin = min_zeta;
max_zeta = max(max(zeta)); p.ColorMax = max_zeta;
% convert seconds since base date into hours since base date
% base_date shoule be set as the cold start date in fort.15 file
time = time/3600; % hour time series since the start time
Base_date = datetime(p.Base_date,'InputFormat','yyyy-MM-dd HH:mm:ss');

nTimes = length(time);

if p.IterStop==-1
p.IterStop=nTimes;
end

%% read into JTWC best track data
tc = readtable(p.BestTrack);
tc_lon = table2cell(tc(:,8));
tc_lat = table2cell(tc(:,7));
tc_lon_lat = zeros(56,2);
for i = 1:length(tc_lon)
tc_lon_t = str2double(cell2mat(extractBefore(tc_lon(i,1),"E")))/10;
tc_lat_t = str2double(cell2mat(extractBefore(tc_lat(i,1),"N")))/10;
tc_lon_lat(i,1) = tc_lon_t;
tc_lon_lat(i,2) = tc_lat_t;
end

%% we use VideoWriter to make the Animation
Animation = VideoWriter([p.VideoWriteDir,'video']);
Animation.FrameRate = 5;
open(Animation)

disp("############")
disp("Start to Generate Image for Each Time Step")
disp("############")
disp("Image is Generating for Date of:")

%% plot image for each timestep and read as one frame
figure
set(gcf,'unit','normalized','position',p.FigurePosition)
set(gca,'FontSize',p.FontSize);
m_proj(p.Projection,'lon',[min(x)-0.25 max(x)+0.25],'lat',[min(y)-0.25 max(y)+0.25]);
m_grid('linestyle','none','box','fancy','tickdir','in');%
title(p.Title,FontSize=p.FontSize+4);
caxis([p.ColorMin p.ColorMax]);
colormap(p.ColorMap);
cb = colorbar;
cb.XLabel.String = 'm'; % {'String','Rotation','Position'}
% loop for each time step
for idx = p.IterStart:p.IterStride:p.IterStop
current_date = Base_date + time(idx,:)/24;
disp(current_date);
current_ele = zeta(:,idx);
m_trisurf(element,x,y,current_ele);
m_line(tc_lon_lat(:,1),tc_lon_lat(:,2),'linewi',2,'color','k');
%m_text(92.5,21.5,datestr(current_date),FontSize=18);
frame = getframe(gcf);
writeVideo(Animation,frame);
end
close(Animation);

disp("Animation done!")

end

[Jiangchao3]

This is just a demo version, I will keep working to make it better.

And, Animation64nc Animation6364nc, Animation6374nc is on the road.

When I complete them, I am happy to push them to OceanMesh2D.

If you have some ideas or great suggestions, please let me know.

[krober10nd]

Hi @Jiangchao3 this looks great! Thanks! we'd love to incorporate this and the additional functionality to the main branch.

[Jiangchao3]

Hi @krober10nd, it is my great pleasure to contribute the OM community.

I am making the animation to show water level and wind stress simultaneously,

use m_quiver to plot the wind in u and v direction.
m_quiver(x,y,current_windx,current_windy,1.5,'color','k','linewidth',0.1);

one problem is that the wind stress values are corresponding to mesh nodes, that will show a very dense vector when it is in nearshore area because there is more high resolution node.

Do you have some suggestions about how to deal with the wind stress values, for example using a resample or other method to generate a uniform wind field, let the visualization looks better?

Thanks

[Jiangchao3]

I have one idea to solve the above wind issues, let me have a try

[krober10nd]

It could be useful to interpolate the data to a structured grid for the quivers. The unstructured grid's variable resolution can make it challenging to see the quivers at times.

[Jiangchao3]

agree with you

[Jiangchao3]

Update: add wind vector into the water level animation:

video.mp4

function [x_grid_reshape,y_grid_reshape,wind_intp_reshape] = interpolant_nc74(x,y,wind,n_grid)
% function to interpete the original unstrural wind into structural grid
% Author: Jiangchao Qiu, (MIT/ESSG; email:qiujch24@mit.edu)
% Created: May 14 2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[x_grid,y_grid] = meshgrid(linspace(min(x),max(x),n_grid),linspace(min(y),max(y),n_grid));
wind_interp = griddata(x,y,wind,x_grid,y_grid);
x_grid_reshape = reshape(x_grid,[],1);
y_grid_reshape = reshape(y_grid,[],1);
wind_intp_reshape = reshape(wind_interp,[],1);
end

but there is some edge error exist after doing the interpolation, should drop them

[Jiangchao3]

Drop the wind interpolation out of the mesh domain:

function [x_grid_indom,y_grid_indom,wind_grid_indom] = interpolant_nc74(x,y,wind,n_grid)
% function to interpete the original unstrural wind into structural grid
%
% Author: Jiangchao Qiu, (MIT/ESSG; email:qiujch24@mit.edu)
% Created: May 14 2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[x_grid,y_grid] = meshgrid(linspace(min(x),max(x),n_grid),linspace(min(y),max(y),n_grid));
wind_interp = griddata(x,y,wind,x_grid,y_grid);
x_grid_reshape = reshape(x_grid,[],1);
y_grid_reshape = reshape(y_grid,[],1);
wind_intp_reshape = reshape(wind_interp,[],1);
%% drop values that are not in the mesh area
% ref: https://www.mathworks.com/help/matlab/ref/inpolygon.html
% ref: https://www.mathworks.com/help/matlab/ref/boundary.html
k = boundary(x,y);
x_boundary = x(k);
y_boundary = y(k);
in = inpolygon(x_grid_reshape,y_grid_reshape,x_boundary,y_boundary);
x_grid_indom = x_grid_reshape(in);
y_grid_indom = y_grid_reshape(in);
wind_grid_indom = wind_intp_reshape(in);
end

[Jiangchao3]

Create a PR here: 5f27b31, please review.

the final version used this function is as following:

video.mp4

[Jiangchao3]

And another version: (add some post-process using Google Earth Studio and Adobe After Effect)

https://www.youtube.com/watch?v=muk0FL4E-Aw&t=15s