matlab2010b有关NCEP数据的个人经验与心得

斥鷃

readnc.m (921 Bytes, 下载次数: 307)

2013-5-2 23:41 上传

点击文件名下载附件
经度取0到357.5的readnc，加投影
不过耐心看完整篇应该用不着这个程序了）
如深深所言，论坛里面有关NC的数据材料信息已经有很多了，不过matlab的方面还是略显老气，这里我再补充一下我使用的2010b的心得体会，算是个手把手的初级教程，用2010b自带的netcdf工具。
   PS：其实ArcMap的多维分析工具也可以读取NC数据，而且不需要进行数据还原。甚至可以进行图层-NC数据的转换，不过略显缓慢。有兴趣的朋友可以试试。
有关NCEP数据的基础内容，NCEP数据的下载，在ArcGis10的中文介绍里面也有有关netcdf的介绍。用grads提取NC数据的方法（链接1，链接2），以及matlab老版的NC读取方法和老版NC工具资源包下载，都可以在论坛里面找到。不过至少在2010版以后的matlab已经没有必要下工具箱读取NC了。前阵子又有坛友发了文章，有关2012版matlab读取NC的方法的帖子，其实我也早想发那么篇贴，不过至今才真正开始动手写，我测试了一下我的这个版本好像没有上面帖子里的ncread这个命令，我就下决心写这么篇帖子。好了，these are previous-works. 下面是我的一些总结和写好的matlab程序（适用于有netcdf工具箱的matlab版本，中文注释）。

首先，几个我们必须知道的NC基本知识（部分参考NCEP数据基础内容）：
    1.NC数据的层次、时次信息和内容可从下载网站上了解得到，例如层次1表示为1000hPa，诸如此类。
    2.NC数据的经纬度、时次等信息会作为数据一部分与我们所要的气温、降水等数据一同存入变量数据中；
    3.NC数据是对原始观测或者说探测数据进行过处理的，通过非专业的软件（地学系统下的grads和arcgis等会自动转换数据）读取的NC数据需要通过计算还原原始值，其公式为：
                                            origindata=readdata*scale_factor+add_offset
其中scale_factor和add_offset数值作为NC数据的属性信息存入属性数据中；
    注记：NC数据存储的2.5*2.5数据的纬度的存储是从90至-90的，与我们一般绘图方式相反，经度则使0度经线为0,东经为0~180的范围，西经则从180算起到357.5，ie. 70W等价于（360-70）=290

好，下面开始进入正题（再增加一部分投影的问题）：
根据上面有关NC数据的论述，我们的读取工作分这么三步：
1.读取NC数据的存储；
2.数据转换（部分数据诸如位势高度，数据不需转换，即scale_factor=1，add_offset=0）；
3.投影问题。

下面以介绍关键matlab命令方式来呈现方法（欲知端的，请听help netcdf分解）：
准备工作，查看数据：
fid=netcdf.open(filename,option)  打开NC数据，返回为一个打开文件的文件号，代表文件；
varnumber=netcdf.inq(fid)           查看文件内变量总数；
读取使用的数据：
varname=netcdf.inqVar(fid,vid)    变量在文件中以vid的形式编号，可以一个个变量看它们的名字，来寻找自己需要的变量，例如多层多年平均气温数据中，vid为4的时候是所需数据变量，0为层次，1为纬度，2为经度，3是时次，而在地表气压数据中，vid为3时为所需数据，其余数据排位都向前顺延，没有层次这一维；
ncdata=netcdf.getVar(fid,vid)       知道自己要找的变量ID号以后，通过getVar命令读取变量（作为变量，经纬度和时次也可以用该命令读取）；
数据转换：
scale_factor=netcdf.getAtt(fid,vid,'scale_factor')   属性数据至今鄙人也不是很清楚具体有哪些，反正scale_factor和add_offset够用了，可以用这条命令读取这两个数据，其中vid是使用变量的id，不要用经纬度的id哦~
origindata=readdata*scale_factor+add_offset      转换成原始数据，数据计算结束。

投影问题：
因为直接读取的NC文件投影方式也直接写在文件里面，因为绘图的关系有必要进行提取。
netcdf.getVar(fid,vid)                                            读取经纬度数据，一般vid=1时是纬度，vid=2时是经度，其中经度有需要的话可以重新进行处理，将>180的数据用公式（data-360）处理为西经的模式，一般画中国地区的话不是很需要这样的转换。
netcdf.close(fid)                                                    最后养成良好习惯，关闭文件。
注记：1.转换那步在计算的时候可能会出现数据类型不匹配的问题，可能需要用到double或者single函数转换数据后方能运算；
           2.要是计算机内存有限啥的（just like my computer~），最后一步计算可能会涉及内存溢出，所以建议先提取所需要的具体数据以后再进行转换，以减小运算规模。
           3.特别注意直接从NC读取出来的数据一般第一维既行标表示的是经度，第二维是纬度，所以绘图的时候按照我们正常习惯需要对提取的数据进行转置处理。

写在最后：
上气象家园一年了，收获灰常灰常大，十分喜欢这个平台，继续四处推广ing。不过久了感觉应该为论坛做些什么，也不光只是下东西。新年新希望，希望论坛越办越好！希望哪天我也可以像大神们那样知道各种各样的东西！
{:soso__9523217983520424030_1:}

斥鷃

deanchen 发表于 2014-9-15 18:30
fid=netcdf.open(filename,option)  打开NC数据，返回为一个打开文件的文件号，代表文件；
varnumber=netc ...

option是代表参数，一般是"NO_WRITE"这种~help 一下netcdf.open

shawjan

感谢楼主分享个人经验，学习！！！

kongfeng0824

不错。支持一下！我的那个新版本的2012b写的nc读取的帖子。里面的函数ncread

斥鷃

kongfeng0824 发表于 2013-1-18 22:12

                               
登录/注册后可看大图

不错。支持一下！我的那个新版本的2012b写的nc读取的帖子。里面的函数ncread

你的文章我也引了，不过我的2010b里help ncread找不到函数，我不知道是不是版本更新的问题，要是有的话用你的那个应该会更简洁些，不过不知道ncread有没有处理了数据还原的问题。

kongfeng0824

斥鷃 发表于 2013-1-18 23:21

                               
登录/注册后可看大图

你的文章我也引了，不过我的2010b里help ncread找不到函数，我不知道是不是版本更新的问题，要是有的话用 ...

嗯。我的版本是2012b的。数据还原你指的是什么意思？

斥鷃

kongfeng0824 发表于 2013-1-18 23:49

                               
登录/注册后可看大图

嗯。我的版本是2012b的。数据还原你指的是什么意思？

像我帖子里说的那样，一部分NC数据为了方便存储是做过处理，我这个版本的话直接读取得到的数据是有问题的，比如读取的700hPa气温数据会是-30000多，所以要还原。不过要是你做出来的图数量级都没什么问题的话应该是对的。

kongfeng0824

斥鷃 发表于 2013-1-18 23:57

                               
登录/注册后可看大图

像我帖子里说的那样，一部分NC数据为了方便存储是做过处理，我这个版本的话直接读取得到的数据是有问题的 ...

嗯。我现在为止还没发现数量级出现差错。多谢！

kongfeng0824

kongfeng0824 发表于 2013-1-18 23:59

                               
登录/注册后可看大图

嗯。我现在为止还没发现数量级出现差错。多谢！

再请教你一个问题：我的数据的全称是CMORPH_V1.0_RAW_0.25deg-DLY_00Z_20080304。windows重命名时默认为该数据的后缀是.25deg-DLY_00Z_20080304。这类数据我看了一下用fopen、fread打开了。但是打开之后只有一列，几十万行的数据。我要把它处理，画在世界地图之上的。不知道你是否有批处理和数据投影的建议呢？

Aires

上面的第二点，add_offset=0吧

Aires

感谢楼主一直支持着论坛，希望能多多贡献好的经验