Two ways to plot a transection

1. Method 1 from Chenfu

clc;clear all;close all;
day=300;
hour=12;
node=6600;
lw=2;

ncload('superior_0001.nc','lon','lat','nv','h','siglay');
lon=lon-360;
ncload(['superior_temp_',num2str(day,'%.4i'),'.nc'],'temp');
map=load('~/Documents/Work/FVCOM_output/Lake_Superior/boundary_superior.dat');

% degree 0
% line_x=[-88.2 -84.7]';
% line_y=[lat(node) lat(node)]';

% degree 90
% line_x=[lon(node) lon(node)]';
% line_y=[46.5 48.6]';

% degree 45&135
degree=45;
% line_x=[-86.21 -87.82]';%135
line_x=[-87.43 -85.89]';%45
k=tan(deg2rad(degree));
line_y=k*line_x+(lat(node)-k*lon(node));

cv=[0:1:18]';

degree=rad2deg(atan((line_y(2)-line_y(1))/(line_x(2)-line_x(1))));

disx=lon(node)-line_x(1);
disy=lat(node)-line_y(1);
distance_node=sqrt(disx.^2+disy.^2);
x_node=[distance_node distance_node];
y_node=[0 -h(node)];

figure(123);hold on;box on;
plot_fvcom_obc(map,[210 180 140]/255);
scatter(lon(node),lat(node));

x=lon;
y=lat;
data=squeeze(temp(hour,:,:));

close all
[out_xyz,topo_xy]=fvcom_section_px_ch(x,y,nv,h,siglay,data,line_x,line_y,cv);
text(0.05,-250,['Day ',num2str(day,'%.3i'),'   Degree ',num2str(degree,'%.3i')],'background','w');
colorbar('location','eastoutside');

% plot(x_node,y_node,'--','Color','r','LineWidth',lw)

saveas(gcf,['fig_transection_deg',num2str(degree,'%.3i'),'_day',num2str(day,'%.3i'),'.png'],'png');  

function [dis,out_x,out_y,out_z,topo_xy]=fvcom_section_px_ch(x,y,nv,h,siglay,data,line_x,line_y,cv)
%[out_x,out_y,out_z,topo_xy]=fvcom_section_px_ch(x,y,nv,h,siglay,data,line_x,line_y,cv)
%function[out_xyz,topo_xy]=fvcom_section_px_ch(x,y,nv,h,siglay,data,line_x,line_y,cv)-%ori
%[out_xyz,topo_xy]=fvcom_section(x,y,nv,h,siglay,data,line_x,line_y,cv)
% out_xyz is the section x(horizontal),y(vertical),z(data)
% topo_xy is the topography xy for land paste.
% cv is the contour vector
% plot section of fvcom data : u ,v ,el, t or s
% section defined by line_x,line_y
% siglay must be layer x node, nv must be 3 x cell, data must be layer x
% node or layer x cell
h=h(:);
xc=mean(x(nv));
yc=mean(y(nv));
%[inter_x,inter_y]=linextriangles(line_x,line_y,nv,x,y); %-ori
[inter_x,inter_y]=linextriangles_ch(line_x,line_y,nv,x,y); %-ch modified to speedup
%inter_x=linspace(line_x(1), line_x(2),100);
%inter_y=linspace(line_y(1), line_y(2),100);

tmp_h=griddata(x,y,h,inter_x,inter_y);
%control section direction
if tmp_h(1)>tmp_h(end)
    inter_x=flipud(inter_x);
    inter_y=flipud(inter_y);
else
    %keep original inter_x,inter_y
end
%end control
[lay,node]=size(siglay);siglay(1,:)=0; 
h_siglay=siglay.*repmat(h',lay,1);
inter_hh=griddata(x,y,h,inter_x,inter_y);
clear inter_h
for i=1:lay
    inter_h(:,i)=griddata(x,y,h_siglay(i,:),inter_x,inter_y);
end
inter_h=inter_h';
%data=h_siglay;
clear inter_data
for i=1:lay
    i
if (length(data)==length(x))
    inter_data(:,i)=griddata(x,y,data(i,:),inter_x,inter_y);
end
if (length(data)==length(xc))
    inter_data(:,i)=griddata(xc,yc,data(i,:),inter_x,inter_y);
end
end
inter_data=inter_data';
disx=inter_x-inter_x(1);
disy=inter_y-inter_y(1);
dis=sqrt(disx.^2+disy.^2); 
dis2=dis;
dis=repmat(dis',lay,1);
dis1=reshape(dis,[],1);
inter_h1=reshape(inter_h,[],1);
inter_data1=reshape(inter_data,[],1);
%[C,H]=scattercontourf(dis1,inter_h1,inter_data1,cv,0);%--chenfu turn off for test
out_xyz=[dis1,inter_h1,inter_data1];
%clabel(C,H,'FontSize',15,'Color','r','Rotation',0)
hold on
%--chenfu turn off for test
%fill([dis2 ;dis2(end:-1:1)],[-inter_hh;-max(inter_hh)*ones(length(inter_hh),1)],'k')
%topo_xy=[[dis2 ;dis2(end:-1:1)],[-inter_hh;-max(inter_hh)*ones(length(inter_hh),1)]];


%figure(33);clf;hold on;
%contourf(flipud(inter_data'));
contourf(dis,inter_h,inter_data,[cv],'linestyle','none');
caxis([cv(1) cv(end)]);

fyy9=[inter_h(end,:) min(inter_h(:)) min(inter_h(:)) inter_h(end,1)];
fxx9=[dis(1,:) dis(1,end) dis(1,1) dis(1,1)];
fill(fxx9,fyy9,'k');

out_x=[dis];
out_y=[inter_h];
out_z=[inter_data];
topo_xy=[fxx9' fyy9'];

2. Method 2 from Chuyan

clc;clear all;close all;

hour=12;
day=250;
interp_layers=[0:10:300]';
line_x=[-88.2 -84.7]';
line_y=[47.1245 47.1245]';
nn=100;
radius=0.5;
minum=1e-5;
npow=2;
nquad=2;
cv=[0:1:18]';

ncload('superior_0001.nc','h','siglay');
ncload(['superior_temp_',num2str(day,'%.4i'),'.nc'],'temp','lon','lat');
lon=lon-360;
[m,n]=size(siglay);

%add surface temp
temp=[temp(:,1,:) temp];
surf_layer=zeros(1,n);
siglay=[surf_layer;siglay];
[m,n]=size(siglay);

h_repmat=repmat(h',m,1);

h_siglay=-h_repmat.*siglay;

%interp temp on every sigma layer

for i=1:n

    X=h_siglay(:,i);
    V=squeeze(temp(:,:,i));
    Xq=interp_layers;
    Vq = interp1(X,V',Xq,'spline','extrap')  ;
    temp_interp_layers (:,:,i) = Vq';
end

%define the position of the transection
inter_x=linspace(line_x(1), line_x(2),nn);
inter_y=linspace(line_y(1), line_y(2),nn);

for i=1:length(interp_layers)

    i
    temp_time=squeeze(temp_interp_layers(hour,i,:));
    [z_obs]=interp_scatter(lon,lat,temp_time,inter_x,inter_y,radius,minum,npow,nquad);
    temp_layer(i,:)=z_obs';

end

disx=inter_x-inter_x(1);
disy=inter_y-inter_y(1);
dis=sqrt(disx.^2+disy.^2); 

%plot
figure(111);hold on;box on;
contourf(dis,-interp_layers,temp_layer,[cv],'linestyle','none');
caxis([cv(1) cv(end)]);

[h_obs]=interp_scatter(lon,lat,h,inter_x,inter_y,radius,minum,npow,nquad);
inter_h=-h_obs';
fyy9=[inter_h(end,:) min(inter_h(:)) min(inter_h(:)) inter_h(end,1)];
fxx9=[dis(1,:) dis(1,end) dis(1,1) dis(1,1)];
fill(fxx9,fyy9,'k');

saveas(gcf,['fig_transection_by_chuyan.png'],'png');