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

652554221

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

参考了一个网上大神的视频。。。只剩下载的视频，找不到出处了，那也得云感谢一下

可以直接看附件，一个文档，先看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.nc是1948年1月-2020年1月

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

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

！！！本质上还是格点法，因为domain度数之后返回的lonk，latk都是点数。。。。

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]);%airbj与bjair48完全一致

%air是三维数据，第一维经度，第二维纬度，第三维是时间。1，12是选取的时间，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开始，也就是四个变量lat、lon、time、air的排序是0、1、2、3。内部数据也是从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.nc是1948年1月-2020年1月

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-864，73个纬度是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_factor和add_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);  %得出时间范围：1948年1月-2015年3月份 （共807个月）；50年1月是第25个数据

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

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

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

上面都是三维数据lon、lat、time，加一个要素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

一大碗年糕

顶一个，干货满满哈哈哈哈

Adiligirl

感谢分享·！！！！！！！！

平凡之路

刚好需要，谢谢

AILSAxl

谢谢楼主！

puck66

刚好需要，谢谢分享

dacheng

宝藏！！真不错，学习了！！！