科研小技巧
从事水文模型或者气候变化相关研究的同学常常会遇到以“.nc”为扩展名的数据文件，不管是前沿文献中公开的各类再分析气象数据还是CMIP6提供的各种全球气候模式输出结果，基本上都会选择以这一格式存储，这就是所谓的“网格通用数据格式”，即NetCDF（以下简称nc）。
图1 CMIP6官网上的模式输出结果，都是以“.nc”为拓展名存储的
nc提供了一个针对高维数据的标准存储格式，被大气、海洋、地球物理等领域广泛应用，但对于接触大气科学较少的诸位水文学家来说怎么查询经纬度，第一次见到nc文件往往不知所措，因为它往往数据量很大，又不能通过常见的通用软件直接打开。高维的数据又比较抽象，打开都做不到，更不要说利用里面的资料了。因此，这里提供一个关于nc的简单介绍，主要是让诸君快速了解nc的基本结构，并且能够读取公开nc资料中的数据，避免诸君了不起的研究被一个小文件给绊倒。
第一步：看看里面有什么？
图2 Panoply从netCDF、HDF、GRIB等数据集绘制地理参考阵列和其他阵列
( )
NASA为查看nc提供了一个免费且高效的工具-Panoply（图2）。其下载和安装详见：
打开Panoply，通过File-Open打开nc文件。以CMIP6的CanESM5模式历史期逐日相对湿度为例，打开后如图3所示。
图3 在Panoply中打开的nc文件
（示例文件为“hurs_day_CanESM5_esm-hist_r1i1p1f1_gn_18500101-20141231.nc”）
图3中左边方框显示本文件中存储的各类数据的列表，包括：
（1）空间坐标：海拔-height、经度-lon、纬度-lat；
（2）时间坐标：time；
（3）变量数据：相对湿度-hurs。
从类型上看，时间和空间坐标都是一维向量（1D）；变量则是一种独特的类型，地理2D（Geo2D）。变量的维度由时间决定，如果时间上只有一点，那么变量数据就是一张地理上2维的地图。
当我们在左框列表中选中某项数据或者数据标题“hurs_day_CanESM5…”时，右框显示数据作者记录在其中的与该项数据有关的各类信息。如果数据作者是个很认真的人，他就会把这部分写得事无巨细，包括数据变量的内容定义“commment”、标准名“standard_name”和曾用名“original_name”、数据的单位“unit”，对于时间，还包括日历采用的格式“calendar”、开始时间、单/双精度，以及数据处理方法、作者名字、处理时间、发布时间、版本等。这些啰里啰唆的信息大部分一般都不会被用到，但当诸君想要充分了解一项不知道从哪里下载来的资料时，就会很感激原作者的细心。
通过Panoly查看，很容易理解nc的数据结构，即由变量矩阵（2维或3维）、时间空间坐标向量以及它们的描述信息组成的一连串数据集。如果诸君熟悉R语言，会发现这个结构跟R语言中的list很相似，都是一串不同类型数据的集合。
除了能打开查看内部信息外，Panoply最重要的功能是对变量矩阵进行可视化。以图3中的相对湿度数据为例，在左框列表中选中“hurs”→点击“Creats Plot”即可逐时间帧绘制相对湿度的空间分布图。
图4 在Panoply中可视化某一天的变量数据
图4即为示例nc文件中的第一帧，也就是hurs矩阵中第一个切片（z=1）数据可视化的结果，它对应的是CanESM5模式对1850年1月1日全球相对湿度的模拟结果。
如果nc的空间坐标是以标准的地理坐标（即经纬度）形式写成的，那么可以就像图4一样通过生成“Longitude-Latitude plot”将你的nc投射到世界地图上；否则需要指定X-Y轴坐标分别是哪个向量，并且只能生成“X-Y plot”。更改“Time：”后面的数字，还可以查看任意帧（即不同时间）的数据图片。这些功能能够帮助我们直观的了解一个nc，并且快速的检查它是否存在比较大的差错。
第二步：读取它！
当然，Panoply只是一个可视化工具，不具备把数据从nc中读取出来并写成我们可用的其他形式的功能，读取并处理nc文件一般需要依赖编程来实现。诸君看到编程先不要紧张。在比较早的时候，读取nc的功能没有得到集成，前辈们只能依靠Fortran等上古语言反复造轮子来处理nc。而在今天，成熟的解释性语言中，基本都已经配置了专门读取nc的工具，例如：R语言的“ncdf4”包，Python的“netCDF4”库，Matlab则内置了“ncread”等一系列处理nc的代码（图5）。这些工具背后的原理不作展开怎么查询经纬度，使用起来大同小异，这里只介绍Matlab中的“ncread”。
图5 Matlab中用于处理nc文件的函数
nc的结构和list类似，实际读取过程也和list类似，需要指定待读取的文件名和变量名。Matlab帮助中心中对于ncread的说明如图6，可见，ncread实际上有三种用法，但一般用第一种也就足够处理绝大部分nc了。第一种用法只需要提供两个参数：source（即想要读取的nc的存储路径+文件名）和varname（即要读取的变量在nc中对应的名字）。
图6 Matlab帮助中心对ncread的说明
图7提供了一个示例，将该文件中的相对湿度“hurs”读出并赋值给Hurs，将经度“lon”读出并赋值给Lon，将纬度“lat”读出并赋值给Lat。因为我们之前通过Panoply已经知道了该文件中各变量的名字，tic和toc则用来计算执行程序需要的时间。
图7 一个读取nc中各变量的Matlab代码示例
图8则显示了上述代码运行后的结果，显然可以得到一个名为Hurs的128×64×60225的三维矩阵、一个名为Lat的64行的向量以及一个名为Lon的128行的向量。需要注意的是，Hurs的行数（128）和列数（64）分别与Lon和Lat的行数一致，说明Hurs的行和列分别对应着经度和纬度，60255则对应着1850年至2014年的天数（按一年365天计算，不考虑闰年，这涉及到CMIP6模式采用的日期计算方法，不赘述）。
图8 示例代码运行后得到的各变量
这样读取得到的变量数据就被以矩阵的形式存储起来，诸位就可以用处理一般矩阵的方式处理它们，比如逐行或者逐列求平均、矩阵乘法、求特征值、检索某一行某一列的数据等等，又或者进一步进行分析，视诸君研究需求所定。
需要提醒的是，上述代码运行一次大约需要20s（具体耗时视读取的nc大小和计算机性能而定），其中比较花时间的是打开nc文件的步骤。这里不做进一步展开，只是提醒诸君，在计算机内存允许的情况下，最好一个文件只读取一次，把读到的信息装起来备用，避免重复多次的读取操作，否则程序运行的时间会很长，白白消耗诸君宝贵的光阴。
补充：检索的思路
本来本篇文章到了第二步已经结束了，但还请容许我就如何从nc中检索特定点的数据多啰嗦几句。
因为nc文件中的变量数据往往是三维的，对于习惯了二维数据的诸位恐怕不那么直观，但只要抓住上文提到的Hurs的行和列分别与经度和纬度（当然，有些nc里会记为X和Y，道理一样）一一对应的关系，即使想不明白“三维数据是个什么名堂？”，也不妨碍诸位进一步工作。
具体来说，Lat中的第一个元素，即对应着Hurs中第一列的经度，Lon中的第一个元素则对应着Hurs中第一行的纬度，其他的也是依次类推。那么当知道一个点的经纬度，如何获得这个点的变量的时间序列呢？聪明的诸君想必很快就想明白了，只要查询该点经纬度分别位于Lat和Lon中的位次，然后依此检索变量中的对应列和行即可。这些操作在不同平台上实现的形式不一，但根本的思路是一致。
当然了，还有一些细枝末节的问题。比如诸位要研究的点的坐标在nc的坐标序列中无法查到，这时候就需要选取研究点临近的数据进行插值，又或者地理坐标和投影坐标的转换问题等等，这些就需要诸君开动脑筋，随机应变了。
希望这篇啰里啰唆的小tips能够为诸位的研究提供一点点微小的帮助，祝愿诸君研究顺利！
——王康铭
供稿| 王康铭
编辑|谈幸燕
校对| 邓皓东
审核| 蒋云钟李玓瓅
出品| 水资源所新媒体工作室
———END———
限 时 特 惠： 本站每日持续更新海量各大内部创业教程，一年会员只需98元，全站资源免费下载 点击网站首页每天更新
站 长 微 信： aiwo51889
