zoukankan      html  css  js  c++  java
  • 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、湿度

      

      

      

      

      

      

      

      

      

     

  • 相关阅读:
    搭建hexo个人博客
    Scanner类使用close()方法问题
    记录一次Ubuntu基础配置和美化
    Linux更换默认Shell
    python-成员修饰符
    ysoserial项目之URLDNS利用分析
    Apereo Cas4.x 反序列化漏洞复现之复现分析与利用
    JAVA反序列化漏洞之调试环境搭建(含ysoserial项目)
    虚拟机window7忘记密码,如何重置?
    多种类型SQL注入
  • 原文地址:https://www.cnblogs.com/ruo-li-suo-yi/p/7762721.html
Copyright © 2011-2022 走看看