各类应用示例

字符串处理

当前路径下包含文件名以'fnl_20110705'开头的4个文件：

 fnl_20110705_00_00_c
 fnl_20110705_06_00_c
 fnl_20110705_12_00_c
 fnl_20110705_18_00_c

fn=dir('fnl_20110705*');
则fn.name中保存了每个文件名字符串

将前两个文件名与第三个文件名前8个字符连接起来：

fns=[fn(1:2).name, fn(3).name(1:8)];
fns=strcat(fn(1:2).name, fn(3).name(1:8)); %与上一行结果相同

垂直方向连接成一个字符矩阵：

fns=char({fn(1:2).name,fn(3).name(1:8)});
fns=strvcat(fn(1:2).name,fn(3).name(1:8));

函数cellfun的使用技巧：

fns={fn.name}：（为包含四个文件名字符串的cell数组）
读出这些文件名中包含的时间信息，保存为数值矩阵：
方法一：
tim=zeros(numel(fns),4)
for i=1:numel(fns)
    tim(i,:)=sscanf(fns{i}, 'fnl_%4d%2d%2d_%2d');
end
方法二：
tim=cellfun(@(x)sscanf(x, 'fnl_%4d%2d%2d_%2d')',fns,'UniformOutput',false)；
tim=cell2mat(tim)；

regexp的使用：

将文件名中的小时数读出来：
tim=regexp(fns,‘\d*(?=_00_c)','match')
%tim=regexp(fns,'[0-9]*(?=_00_c)','match')
tim=[tim{:}];
tim =  (cell类型）
    '00'    '06'    '12'    '18’
将文件名中的’00_c’替换为’30’：
fns2=regexprep(fns,'\d*_c','30')

日期和时间函数

计算1950年1月到2012年12月每月1号相对于1950年1月1日的天数，并将日期转换为英文日期形式：

yr=repmat([1950:2012],[12 1]);
mon=repmat([1:12]', [1 size(yr,2)]);
dy=ones(size(yr));
mydate=[yr(:) mon(:) dy(:)];

daten=datenum(mydate)-datenum([1950 1 1]);

sdate=datestr(datenum(mydate));
% sdate=datestr([mydate zeros(mydate)]);  %与上行结果相同
sdate2=datestr(datenum(mydate), 'mmm dd, yyyy');

如果要使用无闰年的相应计算方法，则可使用MetToolbox中的相应函数：

datenum365, datevec365

各类气象数据文件的读写

NetCDF文件

NetCDF文件的读写方法较多，可以直接利用Matlab的内置函数ncread, ncwrite, nccreate等high-level函数来实现，也可利用netcdf.open, netcdf.inqVar, netcdf.getVar等low-level函数来完成。利用high-level相关函数来读写nc文件的例子代码:
nctest_inner.m

利用MetToolbox中的NetCDF Toolboxc对netcdf文件进行各种操作，更具特色，往往实现起来更为便捷。利用其读取变量的方法：
    nc=netcdf(‘data.nc’);
    var=nc{‘var’}(:);
    close(nc);

读取并创建NETCDF4压缩文件的方法，具体代码见：
nctest_mettb.m

GRIB1/2

利用MetToolbox中的nctoolbox-1.1.1对grib文件中的变量进行读取：
fn='xxx.grib2';
vname='var';
ds=ncgeodataset(fn);
vv=ds.geovariable(vname);
att=cell2struct(vv.attributes(:,2),vv.attributes(:,1),1);

tim=vv.gettimedata();
disp(datestr(tim,'yyyy-mm-dd HH:MM:SS'))
lon=vv.getlondata();
lat=vv.getlatdata();
zz=double(vv.data(:));

Fortran unformatted 顺序存取文件读取

写文件的F90代码如下：
real:: aa(10,20)
aa=reshape((/((i,j=1,10),i=1,20)/),(/10,20/))
open(12,file=‘test.dat’,form=‘unformatted’) ！顺序存取
do i=1,10
  write(12) aa(i,:)
enddo
close(12)

生成的数据文件test.dat大小为880个字节，而10×20×4=800 bytes，多出来的80个字节从何而来？
Fortran在写顺序存取的unformatted文件时，对于real(4)数组，每个write语句在数组前后各有一个4字节的标志位。

MATLAB读取代码如下：
方法一：
fp=fopen('test.dat', 'r');
aa=fread(fp,[22 10], 'float32');
fclose(fp);
aa=aa(2:end-1,:);

方法二：
fp=fopen('test.dat', 'r');
fseek(fp, 4, 'bof');
aa=fread(fp, inf, '20*float32',8);
fclose(fp);
aa=reshape(aa,[20 10]);

GRADS地图文件读写

利用MetToobox中的函数gsmapread和gsmapwrite可实现对GRADS地图文件的读取和写入操作
将grads自带的地图数据文件hires中的地理坐标数据读出来：
[long lat]=gsmapread('hires');

将中国三大河流（长江、黄河和珠江）数据制作成grads的地图数据文件criver：
load griver
[long2 lat2]=polyjoin({s(1:2).long},{s(1:2).lat});
[long lat]=polysplit(long2,lat2);
gsmapwrite('criver',long,lat);

图形绘制及其对象的句柄控制

通过图形和坐标轴句柄对图形和坐标轴位置与大小的确定

对于批量同一类型的绘图来说，可以针对一幅图，手动控制将图形及坐标轴调到合适后，用如下命令获得其位置信息，再将其用于脚本中的设置：
get(gcf,'pos')
ans =        1286         984         742         267
get(gca,'pos')
ans =
        0.1300    0.1100    0.7750    0.8150

设置：
set(gcf,'pos',[1286         984         742         267])
set(gca,'pos',[0.1300    0.1100    0.7750    0.8150])

通过get或findobj来获取子对象句柄，从而对图形对象进行深层操控

有一组10年逐月的气象数据：
 Name       Size            Bytes  Class
 zz        12x10              960  double
现需要画出10年平均的季节变化及每个月10年均值的置信区间样本数n=10，设置信度为95%（alpha=0.05):

Y=mean(zz,2);
E=tinv(1-alpha/2,n-1)*std(zz,0,2)/sqrt(n);
% 则置信区间即为：[Y-E Y+E]
figure
h=errorbar(Y,E,'b','linew',2,'linestyle','-');
set(gca,'fonts',18)

hl=get(h,'child');
set(hl(1),'color','m','marker','o','MarkerfaceColor','y','markers',8)
set(hl(2),'color','r','linew',5)
set(gcf,'color','w')

需要注意的是：对于当前图形，其句柄在利用绘图函数时即可以显式返回，从而通过其能够对图形进行操控。某一图形的子对象，可以通过get(h,'child'）或findobj来获取。但对于不同版本的matlab，由于其不断改进，有些图形对象被严密的封装在相关类之中，因此并不能将其中的某些子对象分离开来，这时通过子对象句柄对图形的控制就会受到极大限制。

填充等值线与风矢量叠加（风矢量图例的绘制）

绘图部分关键代码：
cv=-2.5:0.5:2.5;
figure('color','w'), hold on
[C h hb]=freecontourf(x,y,w,cv);
set(h,'edgecolor','none');
set(findobj(hb,'type','text'),'fonts',18)
h=quiver(x,y,v,w,1.5,'k','linew',1,'maxheadsize',0.03);
vsize=10;
[har ht arsize]=quiverlegend(h,vsize);

绘制风场图例函数quiverlegend中通过句柄来操控图形的相关代码：
hl=findobj(hq,'type','line');
xv=get(hl(1),'XData'); yv=get(hl(1),'YData');
uu=get(hq,'UData'); vv=get(hq,'VData');
har=annotation('arrow',[xt0-arsize xt0],[yt0 yt0],'headstyle','vback3','headlength',5,'headw',4);

quiverlegend的具体代码可以在MetToolbox的m文件中查看。

区域屏蔽（mask）的实现方法

通过inpolygon函数实现：

[lo la]=meshgrid(70:0.5:140,15:0.5:55);
[clong clat]=chngon('g');
[in on]=inpolygon(lo,la,clong,clat);
in=in|on;

figure,hold on
plot(lo,la,'b.','markers',8)
plot(lo(in),la(in),‘r.’,‘markers’,8)
z(~in)=nan;

whos z in  
  Name       Size             Bytes  Class      Attributes
  in        81x141            11421  logical             
  z         81x141            91368  double

这种方法在边缘会产生锯齿空隙。
Grads的mask功能与此类似，但它需要提前准备好mask文件。

 

通过构造闭合多边形区域

假定有若干闭合区域：long和lat为cell类型数组，每个元素为一个闭合区域的经纬度数据（比如为中国大陆，海南和台湾）
%[lo0 la0]为顺时针方向
lo0=[70 70 140 140 70]';
la0=[15 55 55 15 15]';
x=lo0; y=la0;
for i=1:length(long)
    %[long{i} lat{i}]为逆时针方向
    x=[x; long{i}; lo0(1)];
    y=[y; lat{i}; la0(1)];
end
patch(x,y,’c’,'edgecolor','none’)

 
对于中国区域的完美mask功能可通过直接调用MetToolbox中的函数chnmask或chnmak2来实现。

外部应用程序调用

• 由MATLAB脚本文件生成系统可执行文件

• MATLAB调用外部应用程序接口编写及程序编译

• 通过系统命令和临时文件对外部应用程序的调用

由MATLAB脚本文件生成系统可执行文件

mcc 编译命令： MCC Invoke MATLAB to C/C++ Compiler.
MCC [-options] fun [fun2 ...]

第一次使用时须设置编译环境：
mbuild –setup

编译m程序：
mcc –m myprog.m

生成可执行文件：
Windows: myprog.exe
Mac: myprog.app

设置好相关环境变量即可在操作系统下运行该可执行程序。

调用外部应用程序接口编写及程序编译

mex 编译命令:
编译器设置：mex –setup
mex  -v –Llibpath –llibname myprog.c

应用实例：扩展griddata插值方法（不规则数据网格化）

MATLAB内置的griddata函数插值方法较少：'linear’,'cubic’,'natural’,'nearest’,'v4’
通过C语言调用IDL中的GRIDDATA函数，此函数包含了地学中常用的各种离散点（站点）空间插值方法，如KRIGING，INVERSE_DISTANCE等。

函数接口编写（idl_griddata.cpp）：
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
{
   ........................................
    idl_griddata(zi,x,y,z,xi,yi,meth,ns,nx,ny);
}

static void idl_griddata(double zi[], double x[], double        y[], double     z[],
                                         double xi[],double yi[],char* meth, int ns,int nx, int ny)
{
                ........................................
                IDL_ExecuteStr(buf);
                IDL_ExecuteStr("z=double(z)");
                IDL_ExecuteStr("z=transpose(z)");
                ........................................
}

完整的idl_griddata.cpp代码

编译命令:
mex -IC:\RSI\IDL60\external idl_griddata.cpp C:\RSI\IDL60\bin\bin.x86\idl32.lib
生成的mex文件：idl_griddata.mexw64 （windows下）

调用格式：
zi=idl_griddata(x,y,z,xi,yi,meth), 其中meth为插值方法
如：
zi=idl_griddata(x,y,z,lo,la, 'kri');

通过系统命令和临时文件对外部应用程序的调用

系统命令调用方法
（1）系统path设置
（2）系统命令调用
          !cmd，eval
          system

应用实例:
对GRADS中cressman客观分析方法的调用：gradsoacr
Syntax:
  zi=gradsoacr(x,y,z,xi,yi)
  zi=gradsoacr(x,y,z,xi,yi,radii)
  % Default:
  % radii=[10 7 4 2 1];
因为要通过读写文件交换数据，所以这种方法的执行效率要低一些。

gradsoacr.m的具体代码可在MetToolbox中查看。

MetToobox中部分常用函数简介

Lamber投影地图的绘制：maplam

  Plot Lambert Projection map
  Syntax:
  hc=maplam
  hc=maplam(param)
  ......
  hc=maplam(param, iscoast, iscountry, ischina, isriver, isfish, zp, lspace, degtype, cnres, nriver, glstyle)

常规地图及南海图标的绘制：mapnorm，ssaxes

ssaxes
  Plot Chinese South Sea islands
  Syntax:
  hss=ssaxes
  ......
  hss=ssaxes(meth,wa,rp,fdata,widax,iscolor)

可绘制纹理的填充等值线函数freecontourf的使用

help freecontourf
  Plot contourf for uneven contour value and plot texture in contourf
  Syntax:
  [C,h,hb]=freecontourf(z)
  ......
  [C,h,hb]=freecontourf(x, y, z, cv, cmap, ifirst, ndelta, iswbg, bardir, texture, issharp)

区域模式中地图投影转换

deg2dis.m      - Convert Lambert Projection from degree to distance               
dis2deg.m      - Convert Lambert Projection from distance to degree               
deg2grd.m      - Convert Lambert Projection from degree to grid                   
grd2deg.m      - Convert Lambert Projection from grid to degree

代码示例：
param=[151 111 105 35 76 56 15 55 105 60000 60000];
maplam(param)
[xi yi]=meshgrid(1:param(1),1:param(2));
[lo la]=grd2deg(xi,yi,param);
figure, hold on
plot(lo,la,'b.')
axis([35 175 -5 70])
gmapnorm(1,true,false,false,false,0,[10 5],'g',2)

 

各类colormap的应用

Matlab内置的常规colormap

另外matlab也有其它的colormap生成函数，如：polcmap，demcmap

通过MetToolbox函数nclcmap对ncl丰富colormap的调用

Produce ncl colormap
Syntax:
cmap=nclcmap %help nclcmap, and display available cmap name
cmap=nclcmap('all') %display available cmap name and draw them
cmap=nclcmap(cmapname)
cmap=nclcmap(cmapname,ncmap)
cmap=nclcmap(cmapname,ncmap,isdisp)
%cmapname: can be a string name or a indexed number of colormap
%ncmap: number of color to be interped; 0, an original length of colormap

           

使用方法示例：

cv=-5:2:5;
cmap=nclcmap(4,numel(cv)+1); %index=4的colormap
或直接利用colormap的名称：
cmap=nclcmap('BlAqGrWh2YeOrReVi22',numel(cv)+1);
freecontourf(peaks,cv,cmap);

其它示例及完整代码

Taylor图

代码（通过调用MetToolbox函数TaylorDiag实现）和测试数据（模式模拟与观测的降水序列数据，mat文件）

ENSO功率谱分析

代码与测试数据（BNU-ESM模拟结果，mat文件）

投影转换

使用mstruct、mfwdtran和minvtran函数
代码1  代码2

geotiff文件读写（内含投影转换）

代码

WRF模拟嵌套区域绘制

代码

站点观测降水空间内插、填充等等线绘制、中国区域外屏蔽、地图及南海图标绘制

完整代码

带纹理的填充等值线的绘制

代码

多个相邻多边形区域的合并

利用MetToolBox函数uniteprv实现东北三省与内蒙四个区域的合并，生成一个闭合多边形区域

代码

不同投影，不同视角、不同光照与材质下，全球地形的绘制

利用Map Toolbox

代码