用MATLAB畫925hPa位勢高度分布圖

最新推荐文章于 2022-11-22 13:09:09 发布
最新推荐文章于 2022-11-22 13:09:09 发布
阅读量2.9k 收藏 8
点赞数 1
分类专栏： ATMOS MATLAB黑科技
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/ssy8stephy/article/details/50819581
版权

用GFSdata數據畫單月平均位勢高度分布圖

function []=draw(y,m,d,file)
%畫單月平均
%GFSdata

    %--載入toolbox--%
    %setup_nctoolbox
    
    %--計算經緯度對應的idx--%
    [lon1,lon2,lat1,lat2]=findpos(60,140,15,50);%%%%%%%%%LINR caxis & axis
    figure(1)
    
    %--保存平均高度的矩陣--%
    data=zeros(lat2-lat1+1,lon2-lon1+1);%%%%%%%%%%lat,lon
    
    %--用於findmax()--%
    %maxx=10000;
    
    %--讀和處理數據文件--%
    for i =1:d
        i
        if i<10
            name0=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s0%d_00_000',file,y,m,i);
            name6=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s0%d_06_000',file,y,m,i);
            name12=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s0%d_12_000',file,y,m,i);
            name18=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s0%d_18_000',file,y,m,i);
        else
            name0=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s%d_00_000',file,y,m,i);
            name6=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s%d_06_000',file,y,m,i);
            name12=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s%d_12_000',file,y,m,i);
            name18=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s%d_18_000',file,y,m,i);
        end
        
        data=cutadd(name0,data,lon1,lon2,lat1,lat2);
        data=cutadd(name6,data,lon1,lon2,lat1,lat2);
        data=cutadd(name12,data,lon1,lon2,lat1,lat2);
        data=cutadd(name18,data,lon1,lon2,lat1,lat2);
        %{
        maxx=findmax(name0,maxx);
        maxx=findmax(name6,maxx);
        maxx=findmax(name12,maxx);
        maxx=findmax(name18,maxx);
        %}
    end
    data=data./(4*d);
    mydraw(data,name0,lon1,lon2,lat1,lat2);
    
    %--設置顯示範圍--%
    axis([60,140,15,50])%lon,lat
    title(sprintf('20%s-%s',y,m));
    xlabel('lon')
    ylabel('lat')
    
    %--設定colorbar範圍--%
    caxis([600,900]);%650,800
    
    %--保存圖片--%
    print(1,'-dbmp',sprintf('C:\\Users\\kong\\KHIMatlab\\925\\picture\\GFS\\pic4\\SMALLuni20%s-%s',y,m));    
end





function [data]=cutadd(name,data,lon1,lon2,lat1,lat2)
    da=ncgeodataset(name);
    iso=da.geovariable('Geopotential_height_isobaric');
    
    %--選取數據範圍--%
    data1=iso.data(1,4,lat1:lat2,lon1:lon2);%%%%%%%%%%lat,lon
    
    %--二維化--%
    data1=squeeze(data1);
    data=data+data1;    
    da=[];
    iso=[];
    data1=[];
end

function []=mydraw(data,name,lon1,lon2,lat1,lat2)
    da=ncgeodataset(name);
    iso=da.geovariable('Geopotential_height_isobaric');
    
    %--讀取經緯度--%
    pos1=iso.grid_interop(1,4,:,:);
    lat=pos1.lat(lat1:lat2);%%%%%%%%%%
    lon=pos1.lon(lon1:lon2);%%%%%%%%%%    
    
    %--讀取地圖--%
    gx=shaperead('C:\Users\kong\KHIMatlab\ChinaMap\province\bou2_4p.shp','usegeocoords',true);    
    [gca,h]=contourf(lon,lat,data);
    
    %--畫等值線--%
    clabel(gca,h);
    
    %--設定colorbar文字--%
    h=colorbar;    
    set(get(h,'Title'),'string','gpm');   
    geoshow(gx,'FaceColor','none','EdgeColor','w','LineWidth',1.5);
    max(max(data))
    min(min(data))
end

function [lon1,lon2,lat2,lat1]=findpos(lon01,lon02,lat01,lat02)
    name='C:\Users\kong\KHIMatlab\925\data\1201\gfs05d_120101_18_000';
    da=ncgeodataset(name);
    iso=da.geovariable('Geopotential_height_isobaric');
    pos=iso.grid_interop(1,4,:,:);
    lon1=find(pos.lon==lon01)
    lon2=find(pos.lon==lon02)
    lat1=find(pos.lat==lat01)
    lat2=find(pos.lat==lat02)
end

function [maxx]=findmax(name,maxx)
    da=ncgeodataset(name);
    iso=da.geovariable('Geopotential_height_isobaric');
    data1=iso.data(1,4,:,:);
    data1=squeeze(data1);
    if min(min(data1))<maxx
        maxx= max(max(data1))
    end
    da=[];
    iso=[];
    data1=[];
end

用JRA25數據畫多年單月平均位勢高度分布圖 

function []=drawncep(m,d)
%畫多年單月平均
%JAR

    %--載入toolbox--%
    %setup_nctoolbox
    
    %--計算經緯度對應的idx--%
    [lon1,lon2,lat1,lat2]=findpos(95,120,10,30);%%%%%%%%%LINR caxis & axis
    figure(1)
    
    %--保存平均高度的矩陣--%
    data=zeros(lat2-lat1+1,lon2-lon1+1);%%%%%%%%%%lat,lon
    
    %--用於findmax()--%
    %maxx=10000;
    
    %--讀和處理數據文件--%
    for yy=6:11
        if yy<10
            y=sprintf('0%d',yy)
        else
            y=sprintf('%d',yy)
        end
        for i =1:d
            i
            if i<10
                name0=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s0%d00',y,m,y,m,i);
                name6=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s0%d06',y,m,y,m,i);
                name12=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s0%d12',y,m,y,m,i);
                name18=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s0%d18',y,m,y,m,i);
            else
                name0=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s%d00',y,m,y,m,i);
                name6=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s%d06',y,m,y,m,i);
                name12=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s%d12',y,m,y,m,i);
                name18=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s%d18',y,m,y,m,i);
            end

            data=cutadd(name0,data,lon1,lon2,lat1,lat2);
            data=cutadd(name6,data,lon1,lon2,lat1,lat2);
            data=cutadd(name12,data,lon1,lon2,lat1,lat2);
            data=cutadd(name18,data,lon1,lon2,lat1,lat2);
            %{
            maxx=findmax(name0,maxx);
            maxx=findmax(name6,maxx);
            maxx=findmax(name12,maxx);
            maxx=findmax(name18,maxx);
            %}
        end
    end
    data=data./(4*6*d);
    mydraw(data,name0,lon1,lon2,lat1,lat2);
    
    %--設置顯示範圍--%
    axis([95,120,10,30])%lon,lat
    title(sprintf('%s',m));
    xlabel('lon')
    ylabel('lat')
    
    %--設定colorbar範圍--%
    caxis([650,850]);%650,800
    
    %--保存圖片--%
    print(1,'-dbmp',sprintf('C:\\Users\\kong\\KHIMatlab\\925\\picture\\JAR\\pic1\\uni%s',m));    
end





function [data]=cutadd(name,data,lon1,lon2,lat1,lat2)
    da=ncgeodataset(name);
    iso=da.geovariable('Geopotential_height_isobaric');
    
    %--選取數據範圍--%
    data1=iso.data(1,22,lat1:lat2,lon1:lon2);%%%%%%%%%%lat,lon
    
    %--二維化--%
    data1=squeeze(data1);
    data=data+data1;    
    da=[];
    iso=[];
    data1=[];
end

function []=mydraw(data,name,lon1,lon2,lat1,lat2)
    da=ncgeodataset(name);
    iso=da.geovariable('Geopotential_height_isobaric');
    
    %--讀取經緯度--%
    pos1=iso.grid_interop(1,22,:,:);
    lat=pos1.lat(lat1:lat2);%%%%%%%%%%
    lon=pos1.lon(lon1:lon2);%%%%%%%%%%    
    
    %--讀取地圖--%
    gx=shaperead('C:\Users\kong\KHIMatlab\ChinaMap\province\bou2_4p.shp','usegeocoords',true);    
    [gca,h]=contourf(lon,lat,data);
    
    %--畫等值線--%
    clabel(gca,h);
    
    %--設定colorbar文字--%
    h=colorbar;    
    set(get(h,'Title'),'string','gpm');   
    
    %--畫中國地圖--%
    geoshow(gx,'FaceColor','none','EdgeColor','w','LineWidth',1.5);
    
    %--畫世界地圖--%
    geoshow('landareas.shp', 'FaceColor','none','EdgeColor','w','LineWidth',1); 
    max(max(data))
    min(min(data))
end

function [lon1,lon2,lat2,lat1]=findpos(lon01,lon02,lat01,lat02)
    name='C:\Users\kong\KHIMatlab\925\JAR\200605\anl_p.2006050100';
    da=ncgeodataset(name);
    iso=da.geovariable('Geopotential_height_isobaric');
    pos=iso.grid_interop(1,22,:,:);
    lon1=find(pos.lon==lon01)
    lon2=find(pos.lon==lon02)
    lat1=find(pos.lat==lat01)
    lat2=find(pos.lat==lat02)
end

function [maxx]=findmax(name,maxx)
    da=ncgeodataset(name);
    iso=da.geovariable('Geopotential_height_isobaric');
    data1=iso.data(1,22,:,:);
    data1=squeeze(data1);
    if min(min(data1))<maxx
        maxx= max(max(data1))
    end
    da=[];
    iso=[];
    data1=[];
end


KalariKong
关注
  • 1
    点赞
  • 踩
  • 8
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
专栏目录
python画位势高度图_气候变化位势高度
weixin_39683734的博客
01-15 4247
1、位势高度的公式表达位于重力场中某点物体的势能(E)与它本身的质量(m)之比称为这点的位势(Φ)，从物理学中知道重力势能是E=mgz。对空间某点的位势应为在重力场中任一高度上，单位质量空气相对于海平面所具有的位能所表征的高度称为重力位势高度，简称位势高度，通常以位势米为单位。若用Z表示位势高度，则：位势米这个单位是在1947年国际气象组织高空探测委员会多伦多(Toronto)会议上确定的。以“位...
位势函数 matlab程序
01-21
matlab 作业程序 位势函数部分 位势函数实现的matlab代码 可以直接使用 如果不能就把括号去掉
matlab位势高度图,北半球500 hPa位势高度场时空频谱特征分析
weixin_31493023的博客
03-21 1327
Abstract:Using NCEP2 reanalysis data from 1979 to 2004, the daily 500-hPa geopotential height in the Northern Hemisphere (20°N－90°N) has been expanded into double Fourier series and the signals in the...
matlab最佳高度控制,用MATLAB畫925hPa位勢高度分布圖
weixin_34809884的博客
03-21 377
用GFSdata數據畫單月平均位勢高度分布圖function []=draw(y,m,d,file)%畫單月平均%GFSdata%--載入toolbox--%%setup_nctoolbox%--計算經緯度對應的idx--%[lon1,lon2,lat1,lat2]=findpos(60,140,15,50);%%%%%%%%%LINR caxis & axisfigure(1)%--保存...
python画位势高度图_位势高度计算中气压-高度公式的简化及其误差
weixin_31180441的博客
01-15 2264
位势高度计算中气压高度公式的简化及其误差崔喜爱1顾浩2曹云昌1【摘要】摘要：对高空气象学领域而言,位势高度是重要的一个参量，为天气学和气候学业务应用所提供的高度一般是位势高度。位势高度通常采用探空所实时测量的温度、气压和湿度来计算，计算公式为气压-高度公式(简称压高公式)。本文针对压高公式进行了推导分析，并指出公式中简化的由来以及定量计算了简化所引起的误差，分析表明，压高公式中的重力加速度变量若用...
Matlab 气象图绘制实例
08-17
Matlab气象图绘制，有助于了解利用Matlab软件如何进行气象常用图的绘制流程！
func_hpa.rar_HPA _Rapp_ghorbani _matlab HPA_saleh
09-23
标题中的"func_hpa.rar_HPA _Rapp_ghorbani _matlab HPA_saleh"提到了几个关键元素：`HPA`、`Rapp`、`Ghorbani`、`Matlab`以及`HPA_saleh`。这表明压缩包内容可能与一个名为`HPA`（High Power Amplifier）的项目有关...
HPA的matlab代码
05-03
标题中的"HPA的matlab代码"指的是使用MATLAB编程实现的一种特定的算法或方法，很可能是基于高通量数据的聚类分析。HPA通常指的是高通量成像（High-Throughput Phenotyping）或者高功率放大器（High Power Amplifier...
MATLAB天气图
02-16
MATLAB天气图
850hPa_TMPC_Winds_850hPa风的数据绘图分析_metpy_
10-02
在这个示例中，Cartopy被用来在"Contiguous United States (CONUS)"的地理背景下绘制风场和温度分布，提供直观的气象分析视图。 5. **代码实现**：`850hPa_TMPC_Winds.py`这个文件很可能是实现上述功能的主要代码。...
write.ncl
05-11
write.ncl
850hpa
weixin_51886730的博客
07-14 3834
关于850hpa的一些问题 1.为什么看数值预报，要看850hpa场？850hpa是一种什么样的地方？ 850hPa是这样一个地方，它处于海拔1200~1900m之间，平均高度大约在1300~1500m左右，在某种意义上它正好可以代表近地面边界层的顶端，对地面有很强的相关性。看数值要看850hPa是因为它不高也不低，不会高到脱离地面实际，又不会低到被局部小地形地形明显影响，所以在看数值的时候他可以代表近地面附近的情况，当然这个数据是要变化的。 2.什么是高度场？ 高度场有两种，一...
Python用Basemap画位势高度距平图和5880线
mohenxiang的博客
06-09 2944
​ 本文数据为NCEP/NCAR 再分析数据，感觉这个真是比ERA5好用，就是分辨率和内容类型比ERA5差。链接如下：https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html。本文比较初级，就是画了等值线还有一个画出目标区域的线，shp文件网上有很多，想要可以发邮箱在评论区，我把我收藏的都给可以奉上。文中我把一个大文件裁剪成30个小文件，毫无必要，这个可以再优化下。希望能给大家提供帮助。 ​...
matlab画平面风羽图（彩色）
最新发布
miaobo0的博客
11-22 1389
【代码】matlab画平面风羽图（彩色）
matlab画大气环流
miaobo0的博客
11-20 1532
【代码】matlab画大气环流。
MatLab画三维图和等高图
Sun_Hongbin的博客
04-17 4214
a=xlsread('C:\Users\cw\Desktop\每个位置每小时有多少记录.csv'); % a(37:47,1)=4;a(48:59,1)=5;a(60:64,1)=6;a(65:77,1)=7;a(78:92,1)=8; % a(93:99,1)=9;a(100:104,1)=10;a(105:114,1)=11;a(115:125,1)=12;a(23:36,1)=13; % ...
matlab 累积距平,气象统计方法实习报告
weixin_27188815的博客
03-19 2100
气象统计方法实习报告 目录 实习一 求500hPa高度场气候场、距平场和均方差场3 1、资料介绍3 2．要求3 3、实习结果3 1)、FORTRAN源程序3 (2)、grads文件5 (3)、实习结果6 实习二 计算给定数据资料的简单相关系数和自相关系数11 1、资料介绍11 2、要求11 3、实习结果12 (1)、Fortran源程序12 (2)、程序运行结果：14 实习三 分析中国夏季降水线性...
使用MatLab读取grib(grb)和grib2文件的方法
热门推荐
03-30 3万+
在开始读取之前，有必要先了解一下grib文件是什么。以下内容摘自百度百科：GRIB 码是与计算机无关的压缩的二进制编码,主要用来表示数值天气预报的产品资料。现行的GRIB 码版本有GRIB1 和GRIB2 两种格式。GRIB2较之GRIB1具有加大优点而被广泛使用。如：表示多维数据、模块性结构、支持多种压缩方式、IEEE标准浮点表示法等。 目前有两个版本，第一版（GRIB1）和第二版（GRIB
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值