请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 232|回复: 3

[程序设计] nc数据全解读-2 matlab读取方法ncread、netcdf

[复制链接]

新浪微博达人勋

发表于 2020-3-19 19:29:02 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x

nc数据全解读-2  matlab读取方法ncreadnetcdf

参考了一个网上大神的视频。。。只剩下载的视频,找不到出处了,那也得云感谢一下
可以直接看附件,一个文档,先看ncstudy2.m,再看gaoweiduHGT.m

1高级函数读取法ncread

我们在1中用air=ncread('air.mon.mean.nc','air');直接读取air,数据是144*73*865,可以叫格点法高级函数读取(因为简单所以高级。。。。)。通过观察经纬度排列,知道对应经纬度的格点,直接实现了具体时段、具体区域的数据读取。

但有些时候会导致数据过大,占内存严重,可以用domain直接限定ncread读取的经纬度范围,从而缩小数据量,不直接读air了!!!

1.1读取方法同上篇

1.2 domain限定经纬度范围的数据读取,直接读取用的数据。

domain=[115117.5 37.5 40];%限定北京范围东经115-117.5,北纬37.5-40

ncdisp('air.mon.mean.nc');%hgt.mon.mean.nc19481-20201

lat=ncread('air.mon.mean.nc','lat');

lon=ncread('air.mon.mean.nc','lon');

!!!本质上还是格点法,因为domain度数之后返回的lonklatk都是点数。。。。

lonk1=find(lon==domain(1));%选定lon的起始点东经115度对应第47个经度(间隔2.5),

lonk2=find(lon==domain(2));%选定lon的结束点东经117.5度对应第48个,

latk1=find(lat==domain(4));%选定lat的起始点北纬40对应第21个纬度(间隔2.5),

latk2=find(lat==domain(3));%选定lat的结束点北纬37.5度对应第22个纬度,

%time1=find(lat==domain(5));%选定lat的起始点北纬52.5对应第16个纬度(间隔2.5),

%time2=find(lat==domain(6));%选定lat的结束点北纬5度对应第61个纬度,

%开始选取区域

airbj=ncread('air.mon.mean.nc','air',...

    [lonk1 latk1  1],[lonk2-lonk1+1 latk2-latk1+1  12]);%airbjbjair48完全一致

%air是三维数据,第一维经度,第二维纬度,第三维是时间。112是选取的时间,1是开始,第12结束,若12换成inf,则全读865个月,也可以任意选定时间段。

fori=1:1:12;

airbj2=squeeze(airbj(:,:,i));%airbj2是第12个月的北京区域数值,前面11个数值被覆盖了


airbj1948(i,:)=reshape(airbj2,1,4);%将三维变量变为二维,airbj1948和上一篇中beijing1948完全一致。

end

%最后检测missing_value,以防万一

missing_records=ncreadatt('air.mon.mean.nc','air','missing_value');%检测missing_value,以防万一

airbj1948(airbj1948==missing_records)=NaN;%将缺省值补上

size(airbj1948)




2读取的低级函数netcdf

2.1低级函数变量法读取

nc数据全解读-1数据下载,nc简单matlab读取、数据内部结构介绍中3.1可知,air数据是三个维度,四个变量,重点是低级函数中变量是从0开始,也就是四个变量latlontimeair的排序是0123。内部数据也是从0开始,865个月是0-864


程序如下:

clearall

fid=netcdf.open('air.mon.mean.nc','NC_NOWRITE');%读取nc格式数据,'NC_NOWRITE'写保护,不能更改数据。

ncdisp('air.mon.mean.nc');%air.mon.mean.nc19481-20201

latd=netcdf.getVar(fid,0,'double');

lond=netcdf.getVar(fid,1,'double');  

timed=netcdf.getVar(fid,2,'double');

aird=netcdf.getVar(fid,3,'double');%得出温度(共144*73*865个数据)



得到的数据完全同以前的一样,也可以接着用domain限定区域,用法同上,但是值得注意!!!!摘自别处------低级函数使用domain时千万注意。数据是从0开始排序的,865个月是0-86473个纬度是0-72,我的matlab不存在这个问题,也不大懂这一块,还是用ncread好一点。。。。推测可能跟matlab版本有关,旧的存在,新的不存在这个问题。欢迎指正。

低版本matlab低级函数netcdf使用domain程序如下:

domain=[115117.5 37.5 40];%限定北京地区范围E115-117.5,  N37.5-40

lonk3=find(lon==domain(1));%选定lon的起始点东经115度对应第47个经度(间隔2.5),

lonk4=find(lon==domain(2));%选定lon的结束点东经117.5度对应第48个,

latk3=find(lat==domain(4));%选定lat的起始点北纬40对应第21个纬度(间隔2.5),

latk4=find(lat==domain(3));%选定lat的结束点北纬37.5度对应第22个纬度,

airdbj=netcdf.getVar(fid,3,[lonk3-1 latk3-1  0],...

[lonk4-lonk3+1 latk4-latk3+1 12],'double');%


得到的airdbj 与以前的bjair48一致。

file:///C:/Users/Administrator/AppData/Local/Temp/msohtmlclip1/01/clip_image002.jpg

通过下面命令,可能看到最小值最大值不符合常理

size(airdbj)

maxairdbj=max(airdbj(:));

minairdbj=min(airdbj(:));

这是因为scale_factoradd_offset的原因

因为我的scale_factor(系数)是1,所以是正常的

file:///C:/Users/Administrator/AppData/Local/Temp/msohtmlclip1/01/clip_image004.jpg

如果不是1,则

sf=netcdf.getAtt(fid,3,'scale_factor','double');%得到air的系数

addoff=netcdf.getAtt(fid,3,'add_offset','double');%得到air的余数

airdbjzuizhong=airdbj.*sf+addoff;

可以得出最终正确数据


   


2.2,低级函数不详细区分第几变量,简单粗暴法

ncid=netcdf.open('air.mon.mean.nc','NC_NOWRITE');%读取nc格式数据

lonid =netcdf.inqVarID(ncid,'lon');%返回经度信息,0-180是东经,180-360对应着西经的180-0

lon22=netcdf.getVar(ncid,lonid);  %得出全球经度范围在2.5的精度下(0:2.5:357.5),既全球360度划分为144个区域,1-73既东半球(东经0-180度);73-144既西半球(西经180-0度)

latid =netcdf.inqVarID(ncid,'lat');

lat22=netcdf.getVar(ncid,latid);%得出全球纬度范围在2.5的精度下,既全球180度划分为73个区域,1-37既北半球(北纬90-0);37-73既南半球(南纬0-90);

%全球纬度的格式是90-2.5-90,共有73个数据,赤道正好是0;是第37个数据

timeid =netcdf.inqVarID(ncid,'time');%返回时间信息  

time22=netcdf.getVar(ncid,timeid);  %得出时间范围:19481-20153月份 (共807个月);501月是第25个数据

airid =netcdf.inqVarID(ncid,'air');%返回变量850hpa高度的值

air22=netcdf.getVar(ncid,airid);%得出高度(共144*73*807个数据)

3更高维度,四维数据五维变量的读取

上面都是三维数据lonlattime,加一个要素air的读取,有时还会有高度层选取,比如hgt.mon.mean.nc数据有17层,也就是四维数据五维变量

从下面网址选hgt数据

https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.pressure.html  file:///C:/Users/Administrator/AppData/Local/Temp/msohtmlclip1/01/clip_image006.jpg

ncdisp('hgt.mon.mean.nc');

之后发现4维、5变量

Dimensions:

           level = 17

           lat   = 73

           lon   = 144

           time  = 865  (UNLIMITED)

Variables:

    level

           Size:       17x1

           Dimensions:level

          Datatype:   single

           Attributes:

                      units        = 'millibar'

                      long_name    = 'Level'

                      positive     = 'down'

                      GRIB_id      = 100

                      GRIB_name    = 'hPa'

                      actual_range = [1000    10]

                      axis         = 'Z'

    lat  

           Size:       73x1

           Dimensions:lat

          Datatype:   single

           Attributes:

                      units         = 'degrees_north'

                      actual_range  = [90 -90]

                      long_name     = 'Latitude'

                      standard_name = 'latitude'

                      axis          = 'Y'

    lon  

           Size:       144x1

           Dimensions:lon

          Datatype:   single

           Attributes:

                      units         = 'degrees_east'

                      long_name     = 'Longitude'

                      actual_range  = [0         357.5]

                      standard_name = 'longitude'

                      axis          = 'X'

    time

           Size:       865x1

           Dimensions:time

          Datatype:   double

           Attributes:

                      long_name       = 'Time'

                      delta_t         = '0000-01-0000:00:00'

                      avg_period      = '0000-01-0000:00:00'

                      prev_avg_period = '0000-00-01 00:00:00'

                      standard_name   = 'time'

                      axis            = 'T'

                      units           = 'hours since1800-01-01 00:00:0.0'

                      actual_range    = [1297320  1928472]

    hgt  

           Size:       144x73x17x865

           Dimensions:lon,lat,level,time

          Datatype:   single

           Attributes:

                      long_name               = 'Monthlymean geopotential height'

                      valid_range             =[-700  35000]

                      units                   = 'm'

                      add_offset              = 0

                       scale_factor            = 1

                      missing_value           =-9.969209968386869e+36

                      precision               = 0

                      least_significant_digit = 0

                      GRIB_id                 = 7

                      GRIB_name               = 'HGT'

                      var_desc                ='Geopotential height'

                      level_desc              ='Multiple levels'

                       statistic               = 'Mean'

                      parent_stat             = 'Other'

                      dataset                 = 'NCEPReanalysis Derived Products'

                      actual_range            =[-354.458344      32321.0977]

%hgt=ncread('hgt.mon.mean.nc','hgt');直接读取数据量太大144*73*17*865(还得再*8浮点才是数据量),处理较慢

直接domain

lat=ncread('hgt.mon.mean.nc','lat');

lon=ncread('hgt.mon.mean.nc','lon');

lvl=ncread('hgt.mon.mean.nc','level')

time=ncread('hgt.mon.mean.nc','time');

%hgt=ncread('hgt.mon.mean.nc','hgt');

domain=[115 117.5 37.5 40]

lonk1=find(lon==domain(1));%选定lon的起始点东经115度对应第47个经度(间隔2.5),

lonk2=find(lon==domain(2));%选定lon的结束点东经117.5度对应第48个,

latk1=find(lat==domain(4));%选定lat的起始点北纬40对应第21个纬度(间隔2.5),

latk2=find(lat==domain(3));%选定lat的结束点北纬37.5度对应第22个纬度,


%开始选取区域

hgtbj=ncread('hgt.mon.mean.nc','hgt',...

   [lonk1 latk1 3 1],[lonk2-lonk1+1 latk2-latk1+1 1 12]);%第一个中括号应该是开始的经度、纬度、层次、时间,第2个中括号总共要读的经度数、纬度数、总的层次数量、时间数

%[lonk1 latk1 3 1]3代表第三层850hpa,1代表从第一个月开始

%[lonk2-lonk1+1 latk2-latk1+1 1 12],1代表只有850一层,12表示到第12个月结束,也可以换成inf,代表到2020年1月

hgtbj2=squeeze(hgtbj)%4维变3维

若继续变2维,跟以前多次处理一致,看下面


for i=1:1:12;

hgtbj3=squeeze(hgtbj(:,:,i));%hgtbj3是第12个月的北京区域850hap数值,前面11个数值被覆盖了


hgtbj48(i,:)=reshape(hgtbj3,1,4);%将三维变量变为二维,


end

gaoweiduHGT.m

1.22 KB, 下载次数: 2, 下载积分: 金钱 -5

ncstudy2.m

3.08 KB, 下载次数: 3, 下载积分: 金钱 -5

nc数据全解读-2 matlab读取方法ncread、netcdf.docx

290.35 KB, 下载次数: 5, 下载积分: 金钱 -5

密码修改失败请联系qq:937062711

新浪微博达人勋

发表于 2020-3-19 22:34:36 | 显示全部楼层
顶一个,干货满满哈哈哈哈
密码修改失败请联系qq:937062711
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-3-26 11:08:49 | 显示全部楼层
感谢分享·!!!!!!!!
密码修改失败请联系qq:937062711
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 4 天前 | 显示全部楼层
刚好需要,谢谢
密码修改失败请联系qq:937062711
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表