Matlab——m_map指南(4)——实例

1、 全球/地区温度图

(1)读取数据

clear all
setup_nctoolbox %调用工具包
tic %计时
%% 
nc=ncgeodataset('tmpsfc.gdas.199401.grb2');   %读文件
tem_1=nc.variables   %浏览数据类型
%% 
a1=nc.geovariable(tem_1(1));%取得数据类型为Temperature_surface的数据
b1=a1.data(1,:,:); %第一个时间点温度数据
c1=squeeze(b1)-273.16;%删除单一维度,换为摄氏温度
%% 
a2=nc.geovariable(tem_1(2));%取得数据类型为lat的数据,纬度
b2=a2.data(:,1)%提取数据
%% 
a3=nc.geovariable(tem_1(3));%经度
b3=a3.data(:,1)%
%% 
a4=nc.geovariable(tem_1(4));%取得数据类型为time的数据
b4=a4.data(:,1)%
%% 
save A b2 b3 c1
toc

读取的是NCEP/CFSR数据,1994年1月的温度数据。

该数据一共四项。温度,纬度,经度,时间。温度中是3维数据组织形式

时间数据

经度数据

纬度数据

可以看出该数据集的数据组织形式,经度纬度构成世界地图,记录了744个时间点的温度数据。1小时采集一次数据。

温度数据

保存为mat格式数据,为画图做准备。

(2)画图

clear all
load A

[Plg,Plt]=meshgrid(b3',b2');%形成网格

 m_proj('hammer-aitoff','clongitude',-150);%投影模式

m_pcolor(Plg,Plt,c1);
shading flat;
colormap('jet');%颜色选择
hold on;
m_pcolor(Plg-360,Plt,c1);
shading flat; %着色模式
colormap('jet');

m_coast();
m_grid('xaxis','middle');

% h=colorbar('h');
% set(get(h,'title'),'string','1991年1月全球温度');

c=colorbar('southoutside','fontsize',12)
c.Label.String = '1994年1月全球温度';
c.Label.FontSize = 15;

(3)找出最大、最小温度的经纬度

clear all
load C
Max_col=max(c1);%列最大值
Max_row=max(c1,[],2)%行最大值
Max=max(max(c1));
[x1,y1]=find(c1==max(max(c1)));%x 行,y 列
T_1=Plt(x1,y1)%纬度
T_2=Plg(x1,y1)%经度

Min_col=min(c1);%列最大值
Min_row=min(c1,[],2)%行最大值
Min=min(min(c1));
[x,y]=find(c1==min(min(c1)));%x 行,y 列
T_x=Plt(x,y)%纬度
T_y=Plg(x,y)%经度

 

可以看出最热52度,在澳大利亚那块(142.8123E,23.2610S);最冷-62度,在北极圈那块(89.9999E,66.3486N)。

(4)中国(地区)温度图

clear all
load A

LATLIMS=[3 54];
LONLIMS=[72 134];%选定边界范围

[Plg,Plt]=meshgrid(b3',b2');%形成网格

 m_proj('lambert','lon',LONLIMS,'lat',LATLIMS);%投影模式

m_pcolor(Plg,Plt,c1);
shading flat;
colormap('jet');%颜色选择
hold on;
m_pcolor(Plg-360,Plt,c1);
shading flat; %着色模式
colormap('jet');

m_coast();
m_grid('box','fancy','tickdir','in');

% h=colorbar('h');
% set(get(h,'title'),'string','1991年1月全球温度');

c=colorbar('southoutside','fontsize',12)
c.Label.String = '1994年1月中国温度';
c.Label.FontSize = 15;

该方法是读取全球数据,只展示部分

(5)

改进的区域方法,读取该区域的数据数据,

clear all
setup_nctoolbox %调用工具包
tic %计时
%% 
min_lat=115;
max_lat=279;
min_lon=231;
max_lon=430;  %区域经纬度范围,在数据中的位置

%% 
nc=ncgeodataset('tmp2m.gdas.199401.grb2');   %读文件
tem_1=nc.variables   %浏览数据类型
%% 
N1=nc.size(tem_1(1));%读取数据大小,可以看出数据的组织形式
a1=nc.geovariable(tem_1(1));%取得数据类型为Temperature_surface的数据
b1=a1.data(1,1,min_lat:max_lat,min_lon:max_lon); %第一个时间点温度数据
c1=squeeze(b1)-273.16;%删除单一维度,换为摄氏温度
%% 
N2=nc.size(tem_1(2));
a2=nc.geovariable(tem_1(2));%取得数据类型为lat的数据,纬度
b2=a2.data(min_lat:max_lat,1);%提取数据
%% 
N3=nc.size(tem_1(3))
a3=nc.geovariable(tem_1(3));%经度
b3=a3.data(min_lon:max_lon,1);%
%% 
N4=nc.size(tem_1(4))
a4=nc.geovariable(tem_1(4));%取得数据类型为time的数据
b4=a4.data(:,1);%
%% 
N5=nc.size(tem_1(5));%读取数据大小
a5=nc.geovariable(tem_1(5));%取得数据类型为time的数据
b5=a5.data(:,1);%
%% 
save tem b2 b3 c1
toc

clear all
load tem

LATLIMS=[3 54];
LONLIMS=[72 134];%选定边界范围

[Plg,Plt]=meshgrid(b3',b2');%形成网格

 m_proj('lambert','lon',LONLIMS,'lat',LATLIMS);%投影模式

m_pcolor(Plg,Plt,c1);
shading flat;
colormap('jet');%颜色选择
hold on;
m_pcolor(Plg-360,Plt,c1);
shading flat; %着色模式
colormap('jet');

m_coast();
m_grid('box','fancy','tickdir','in');

% h=colorbar('h');
% set(get(h,'title'),'string','1991年1月全球温度');

c=colorbar('southoutside','fontsize',12)
c.Label.String = '1994年1月中国温度';
c.Label.FontSize = 15;

(6)读取方式的改变

clear all
setup_nctoolbox %调用工具包
tic %计时
%% 
min_lat=115;
max_lat=279;
min_lon=231;
max_lon=430;  %区域经纬度范围,在数据中的位置

%% 
nc=ncgeodataset('tmp2m.gdas.199401.grb2');   %读文件
tem_1=nc.variables   %浏览数据类型

%% 
N1=nc.size(tem_1(1));%读取数据大小
b1=nc.data(tem_1(1),[1,1,min_lat,min_lon],[1,1,max_lat,max_lon]);%初始的读取位置,最终的位置
c1=squeeze(b1)-273.16;
%% 
N2=nc.size(tem_1(2));%读取数据大小
b2=nc.data(tem_1(2),[min_lat],[max_lat]);
%% 
N3=nc.size(tem_1(3));%读取数据大小
b3=nc.data(tem_1(3),[min_lon],[max_lon]);
save tem b2 b3 c1
toc

可以看出海面2米的温度组织形式,4维数据,包含了744个时间点,1个位置,576个纬度点,1152个经度点。

读取方式是初始位置用一个数组表示,终止位置用一个数组表示

clear all
load tem

LATLIMS=[3 54];
LONLIMS=[72 134];%选定边界范围

[Plg,Plt]=meshgrid(b3',b2');%形成网格

 m_proj('lambert','lon',LONLIMS,'lat',LATLIMS);%投影模式

m_pcolor(Plg,Plt,c1);
shading flat;
colormap('jet');%颜色选择
hold on;
m_pcolor(Plg-360,Plt,c1);
shading flat; %着色模式
colormap('jet');

m_coast();
m_grid('box','fancy','tickdir','in');

% h=colorbar('h');
% set(get(h,'title'),'string','1991年1月全球温度');

c=colorbar('southoutside','fontsize',12)
c.Label.String = '1994年1月中国温度';
c.Label.FontSize = 15;

(7)多个时间的平均值

clear all
setup_nctoolbox %调用工具包
tic %计时
%%
min_lat=115;
max_lat=279;
min_lon=231;
max_lon=430;  %区域经纬度范围,在数据中的位置
 
%%
nc=ncgeodataset('tmp2m.gdas.199401.grb2');   %读文件
tem_1=nc.variables   %浏览数据类型
N1=nc.size(tem_1(1));%读取数据大小
d=zeros(1,165,200);%预定义最后的数值存放空间
f=0%验证预留数
%%
for n=1:10  %选取10个时间点
    b1=nc.data(tem_1(1),[n,1,min_lat,min_lon],[n,1,max_lat,max_lon]);%初始的读取位置,最终的位置
    c(n,:,:)=squeeze(b1)-273.16;
    d=d+c(n,:,:); %最终结果
end
for n=1:10
    f=f+c(n,1,1);
end
e=squeeze(d);
% save tem 
toc

f值等于e的第一个值,说明计算正确.最后计算平均值即可。 

2、风向

(1)数据读取

clear all
setup_nctoolbox
tic
%% 读取数据文件
wind= ncgeodataset('wnd10m.gdas.199401.grb2');
wind_list = wind.variables;%文件的列表情况
%% 
size_of_u = wind.size(wind_list(1));%u分量的数据尺寸,777小时,1个高度,经纬度数据,4D数据
data_u=wind.geovariable(wind_list(1));%取得数据类型为风速u的数据
u_1=data_u.data(1,1,:,:); %
u_2=squeeze(u_1);
%% 
size_of_v = wind.size(wind_list(2));%v分量的数据尺寸,777小时,1个高度,经纬度数据,4D数据
data_v=wind.geovariable(wind_list(2));%取得数据类型为风速v的数据
v_1=data_v.data(1,1,:,:); %
v_2=squeeze(v_1);
%% 
size_of_h= wind.size(wind_list(5));%v分量的数据尺寸,777小时,1个高度,经纬度数据,4D数据
data_h=wind.geovariable(wind_list(5));%取得数据类型为风速v的数据
v_1=data_h.data(1); %高度10米
%% 
wind_speed=sqrt(u_2.^2+v_2.^2);  %矢量合成
save wind u_2 v_2 wind_speed
toc

(2)展示

clear all
load lon_lat %载入坐标数据
load wind  %载入风速数据

LATLIMS=[3 15];
LONLIMS=[115 134];%选定边界范围

m_proj('lambert','lon',LONLIMS,'lat',LATLIMS);%投影模式

min_lat=115;
max_lat=279;
min_lon=231;
max_lon=430; 
 
m_coast;
m_grid('box','fancy','tickdir','in');%边缘经纬度宽度
hold on;%图层合并
m_quiver(Plg(min_lat:max_lat,min_lon:max_lon),Plt(min_lat:max_lat,min_lon:max_lon),...
    u_2(min_lat:max_lat,min_lon:max_lon),v_2(min_lat:max_lat,min_lon:max_lon)...
    ,7.5,'r');%经度,纬度,向东,向北的速度分量,转化为矢量箭头
xlabel('global surface winds 1994/01');

(3)

clear all
load lon_lat
load wind

m_proj('hammer-aitoff','clongitude',-150);%投影模式

m_quiver(Plg,Plt,u_2,v_2,15,'r');
shading flat;
colormap('jet');%颜色选择
hold on;
m_quiver(Plg-360,Plt,u_2,v_2,15,'r');
shading flat; %着色模式
colormap('jet');

m_coast();
m_grid('xaxis','middle');

xlabel('global surface winds 1994/01');

  

3、气压

(1)一个时间点的气压

clear all
clear all
setup_nctoolbox
tic
%% 读取数据文件
pre= ncgeodataset('pressfc.gdas.199401.grb2');
pre_list = pre.variables;%文件的列表情况
%% 
size_of_pre = pre.size(pre_list(1));%620小时,576纬度,1152经度
pre_1=pre.data(pre_list(1),[1,1,1],[1,576,1152]); %
pre_2=squeeze(pre_1);
%% 
size_of_2 = pre.size(pre_list(2));%v分量的数据尺寸,777小时,1个高度,经纬度数据,4D数据
size_of_3 = pre.size(pre_list(3));
size_of_4 = pre.size(pre_list(4));
size_of_5 = pre.size(pre_list(5));
size_of_8 = pre.size(pre_list(8));
t_4=pre.data(pre_list(4),1,620); 
t_8=pre.data(pre_list(8),1:124);
toc

可以看出时间点交错开了,所以最终处理数据时还要按顺序排列起来。

 

clear all
load lon_lat
load pre
 m_proj('hammer-aitoff','clongitude',-150);%投影模式

m_pcolor(Plg,Plt,pre_2);
shading flat;
colormap('jet');%颜色选择
hold on;
m_pcolor(Plg-360,Plt,pre_2);
shading flat; %着色模式
colormap('jet');

m_coast();
m_grid('xaxis','middle');

% h=colorbar('h');
% set(get(h,'title'),'string','1991年1月全球温度');

c=colorbar('southoutside','fontsize',12)
c.Label.String = '1994年1月全球气压';
c.Label.FontSize = 15;

 (2)

clear all
setup_nctoolbox
tic
%% 读取数据文件
pre= ncgeodataset('pressfc.gdas.199401.grb2');
pre_list = pre.variables;%文件的列表情况
size_of_1 = pre.size(pre_list(1));%620小时,576纬度,1152经度
size_of_2 = pre.size(pre_list(2));%v分量的数据尺寸,777小时,1个高度,经纬度数据,4D数据
size_of_3 = pre.size(pre_list(3));
size_of_4 = pre.size(pre_list(4));
size_of_5 = pre.size(pre_list(5));
size_of_8 = pre.size(pre_list(8));

t_4=pre.data(pre_list(4),1,620); 
t_8=pre.data(pre_list(8),1:124);

pre_sum=zeros(744,576,1152);
%% 
for n=1:744
    if mod(n-1,6)==0  %取余数
        a=floor(n/6);
        n_time=n-a*5;%下标位置
        pre_sum(n,:,:)=pre.data(pre_list(5),[n_time,1,1],[n_time,576,1152]); %
    else
        a=floor((n-1)/6);
        n_time=n-a-1;%下标位置
        pre_sum(n,:,:)=pre.data(pre_list(1),[n_time,1,1],[n_time,576,1152]);
    end  
end  
%% end

pre_data=squeeze(mean(pre_sum))%求平均值

save pre_data_199401 pre_data
toc

综合运用余数和向下取整的方法,每隔五个从另外的数据中取一个数。最后用mean()函数求平均值。

4、湿度

 

  

  

 

  

 

 

  

  

  

 

  

  

  

 

posted on 2017-10-31 17:41  箬笠蓑衣  阅读(17505)  评论(8编辑  收藏  举报