MATLAB画NCEP风场

clc
clear
ncFilePath='G:\NCEP\u.nc';
ncFilePath2='G:\NCEP\v.nc';
% ncFilePath3='G:\NCEP\omega.nc';
ncFilePath4='G:\NCEP\hgt.nc';
ncFilePath5='G:\NCEP\airT.2018.nc';
s=datenum('01/01 08:00:00');e=datenum('04/18 02:00:00');d=21600/3600/24;TT=datestr((s:d:e)')
T=datenum(TT);

%% 显示结构
ncdisp(ncFilePath);%显示nc文件的所有结构，以便大概了解里面的内容
u=ncread(ncFilePath,'uwnd');%显示指定变量的内容，注意一定要是变量variables才可以
v=ncread(ncFilePath2,'vwnd');
% omega=ncread(ncFilePath3,'omega');
hgt=ncread(ncFilePath4,'hgt');
air=ncread(ncFilePath5,'air');
% ncdisp(ncFilePath,'/','min');%简单显示结构以及定义
% ncdisp(ncFilePath,'/','full');%全部显示所有结构和定义信息
lon=ncread(ncFilePath,'lon');%读取经度变量
lat=ncread(ncFilePath,'lat');%读取纬度变量
lon=lon-180;
level=ncread(ncFilePath,'level');%读取纬度变量
[XI,YI]=meshgrid(lon,lat);
u850=u(1:144,1:73,3,1);
v850=v(1:144,1:73,3,1);
air850=air(1:144,1:73,3,1);
%  omega850=omega(1:144,1:73,3,1);
hgt850=hgt(1:144,1:73,3,1);
%  air850=air(1:144,1:73,3,1);
%  pcolor(XI,YI,hgt850');
%  shading interp;
 figure
[C,h]=contour(XI,YI,hgt850'./10,[136:4:156],'LineWidth',1.2,'ShowText','on','Color','k');
% clabel(C,h,'fontsize',16,'color','r') 
hold on
contour(XI,YI,air850'-273,5,'LineWidth',1.2,'ShowText','on','Color','r');
% clabel(C,h,'fontsize',16,'color','r') 
 hold on
 [x,y]=meshgrid(lon,lat)
  u=u850'.*sin(y);
  quiver(x,y,u850',v850')
%   contour(XI,YI,hgt850',9,'k-','LineWidth',1.2,'ShowText','on','Color','k');
  
 figure
 [x,y]=meshgrid(-1:.2:1,-1:.2:1)
 u=exp(x).*sin(y);
  v=exp(y).*sin(x);
  quiver(x,y,u,v)
 
 t=u(1:144,1:73,1,1);
 pcolor(XI,YI,t');
 shading interp;
 contour(XI,YI,t',9,'k-');
 t1=u(1:144,1:73,3,1);
 figure
  pcolor(XI,YI,t');
 shading interp;
 hold on
 contour(XI,YI,t1',9,'k-');
 %%画风场设置
%  重点1：
% m_quiver(x(1:d:end,1:d:end),y(1:d:end,1:d:end),ustr_1_2(1:d:end,1:d:end)'./dd,vstr_2_1(1:d:end,1:d:end)'./dd,0);
% m_quiver(30,-80,500./dd.*cosd(-80),0,0)
% 为了使得风场和标尺的大小一致，在m_quiver中最后一个值设置为“0”，即令scale=0
% quiver(x,y,u,v,scale);
% 
% 重点2：
%  d = 3;dd =10;
%  m_quiver(x(1:d:end,1:d:end),y(1:d:end,1:d:end),ustr_1_2(1:d:end,1:d:end)'./dd,vstr_2_1(1:d:end,1:d:end)'./dd,0);
% 
% 通过d来调整箭头间距，dd来调整箭头长度（数值大小）
% 
% 重点3：
%  ustr_1_2=ustr_1_1.*cosd(y');
%  h=m_quiver(30,-80,500./dd.*cosd(-80),0,0);
% m_map 中用的是球坐标系，所以要在所有的纬向量上面乘以cosd（纬度）  注意：是cosd（纬度） cos（弧度）
% 其中纬度是每个纬向量所对应的纬度
% 在ustr_1_2=ustr_1_1.*cosd(y');里面，乘以cosd（y'） 要仔细检查，可以一步一步对应着去找。
%*********************************************************************
clear all;close all;clc
 % 利用NCEP_Real-time_Marine_Data的风应力数据 
 %（http://www.esrl.noaa.gov/psd/data/gridded/data.nmc.marine.html）
 % 绘制全球风场分布
 
 
 %地形
 m_proj('miller','lat',[-90,90],'lon',[0,360]);
 m_coast('patch',[.6 .6 .6],'edgecolor','none');
 m_grid('box','fancy','linestyle','none');
 hold on
 
 %NCEP mean风场数据
 datapath_1='E:dataNCEPNCEP_Real-time_Marine_DataU-wind_Stressustr.mean.nc';
 ncdata_1=netcdf.open(datapath_1,'NC_NOWRITE');
 
 lat_1=netcdf.getVar(ncdata_1,0);
 lon_1=netcdf.getVar(ncdata_1,1);
 time_1=netcdf.getVar(ncdata_1,2);
 ustr_1=netcdf.getVar(ncdata_1,3);
 ustr_1_1=ustr_1(:,:,200);
 
 datapath_2='E:dataNCEPNCEP_Real-time_Marine_DataV-wind_Stressvstr.mean.nc';
 ncdata_2=netcdf.open(datapath_2,'NC_NOWRITE');
 
 lat_2=netcdf.getVar(ncdata_2,0);
 lon_2=netcdf.getVar(ncdata_2,1);
 time_2=netcdf.getVar(ncdata_2,2);
 vstr_2=netcdf.getVar(ncdata_2,3);
 vstr_2_1=vstr_2(:,:,200);
 
 [x,y]=meshgrid(lon_1,lat_1);
 ustr_1_1(ustr_1_1<-1e+6)=nan;
 ustr_1_2=ustr_1_1.*cosd(y');
 vstr_2_1(vstr_2_1<-1e+6)=nan;
 d = 3;dd =10;
 m_quiver(x(1:d:end,1:d:end),y(1:d:end,1:d:end),ustr_1_2(1:d:end,1:d:end)'./dd,vstr_2_1(1:d:end,1:d:end)'./dd,0);
 hold on
 h=m_quiver(30,-80,500./dd.*cosd(-80),0,0);
 set(h,'color','r','linewidth',1);
 hold on
 m_text(30,-75,'50m^2/s^2','fontsize',12);
 
 print(gcf,'-djpeg','-r300','E:quiver_ncep');