matlab nc数据绘制位密平面图和剖面图

注：这是我老师给的程序，自己稍微改动了下

clear all;close all;clc
%%
% change the directory according to your data
your_dir = '/Users/boopstore3/Desktop/matlab/shiyanchengxu/seven/';
fname    = [ your_dir 'hawaii_Hydrobase_soest.nc' ];
% read data
lon      = ncread( fname , 'longitude' );
lat      = ncread( fname , 'latitude' );
depth    = ncread( fname , 'depth');
temp     = ncread( fname , 'temp' );
salt     = ncread( fname , 'salt' );
%
[lon2d lat2d] = meshgrid(lon,lat);
lon2d   = lon2d';
lat2d   = lat2d';
depth3d = reshape(depth,[1 1 length(depth)]);
depth3d = repmat(depth3d,[length(lon) length(lat) 1]);

pden = zeros(length(lon), length(lat), length(depth));
% calculate potential density using GSW-seawater
% %%%% GSW Oceanographic Toolbox
% Calculates sea pressure, dep<0
for zi = 1:length(depth)
    pres = gsw_p_from_z( -depth3d(:,:,zi), lat2d);
    % Calculates Absolute Salinity from Practical Salinity.
    [SA, in_ocean] = gsw_SA_from_SP(salt(:,:,zi), pres, lon2d, lat2d);
    % Calculates Conservative Temperature from potential temperature
    % [CT]=gsw_CT_from_pt(SA, temp);
    [CT]=gsw_CT_from_t(SA,temp(:,:,zi), pres);
    % Calculates potential density anomaly with reference pressure of 0 dbar
    pden(:,:,zi) = gsw_sigma0(SA, CT);
end


%%    绘制位密的平面图
figure(1)
set( gcf , 'color' , 'w' );
m_proj('miller' , 'lon' , [ lon(1) lon(end) ] ,  'lat' , [ lat(1) lat(end) ] );
% draw data
[ C1 , h1 ] = m_contourf( lon2d , lat2d , pden(:,:,1) );
h1.LevelList = [ 20: 0.5 : 30 ];
%h1.LineStyle = 'none';
caxis( [ 20 30 ] );

colormap( 'jet' );
h_cbar = colorbar;
h_cbar.FontSize = 13;
h_cbar.Ticks = [20 :2 : 40 ];
h_cbar.Label.String = ' ';

% add land coast
m_coast( 'patch' , [ .7 .7 .7 ] );
m_grid('box','fancy')
title( 'surface potential density' , 'fontsize' , 18 , 'fontweight' , 'bold' );

%%  剖面图
clear all;close all;clc
% -------------------------------------------------------------------------
%%
% change the directory according to your data
your_dir = '/Users/boopstore3/Desktop/matlab/shiyanchengxu/seven/';
fname    = [ your_dir 'hawaii_Hydrobase_soest.nc' ];
% read data
lon      = ncread( fname , 'longitude' );
lat      = ncread( fname , 'latitude' );
depth    = ncread( fname , 'depth');
temp     = ncread( fname , 'temp' );
salt     = ncread( fname , 'salt' );
%
[lon2d lat2d] = meshgrid(lon,lat);
lon2d   = lon2d';
lat2d   = lat2d';
depth3d = reshape(depth,[1 1 length(depth)]);
depth3d = repmat(depth3d,[length(lon) length(lat) 1]);

pden = zeros(length(lon), length(lat), length(depth));
% calculate potential density using GSW-seawater
% %%%% GSW Oceanographic Toolbox
% Calculates sea pressure, dep<0
for zi = 1:length(depth)
    pres = gsw_p_from_z( -depth3d(:,:,zi), lat2d);
    % Calculates Absolute Salinity from Practical Salinity.
    [SA, in_ocean] = gsw_SA_from_SP(salt(:,:,zi), pres, lon2d, lat2d);
    % Calculates Conservative Temperature from potential temperature
    % [CT]=gsw_CT_from_pt(SA, temp);
    [CT]=gsw_CT_from_t(SA,temp(:,:,zi), pres);
    % Calculates potential density anomaly with reference pressure of 0 dbar
    pden(:,:,zi) = gsw_sigma0(SA, CT);
end

%%          位密剖面图
[latgrid depgrid] = meshgrid(lat,-depth);
latgrid = latgrid';
depgrid = depgrid';
pden_sec = squeeze( pden(81,:,: ) );
%
figure(1)
set( gcf , 'color' , 'w' );
% picture position
lef0 = 0.1; bot0 = 0.1; wid0 = 0.85; hei0 = 0.65;
ax  = axes('position',[lef0 bot0 wid0 hei0]);
%
[c2 h2] = contourf(latgrid, depgrid, pden_sec);
h2.LevelList = [ 20: 0.5 : 30 ];
%h1.LineStyle = 'none';
caxis( [ 20 30 ] );
ax.YLim = [-1000 0];
% setting colorbar
colormap( 'jet' );
h_cbar = colorbar;
h_cbar.FontSize = 13;
h_cbar.Ticks = [ 20 : 2 : 40 ];
h_cbar.Label.String = ' ';
title( 'potential density profile' , 'fontsize' , 18 , 'fontweight' , 'bold' );