为了更好地对NCL有一个整体的认知，首先，举一个最简单的例子。
画出1948年以来500hPa高度上的平均等位势线。
数据可从以下网址中获取：https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html
Fig.1.hgt.mon.mean.nc数据下载.png
有了数据，知道了要画什么，首先你要做的是：知道数据的存储格式，也就是知道数据长什么样子。
像我们刚刚下载的数据，属于比较正规的数据，其网站上面有一些说明：
Fig.2.数据格式说明.png
Brief Description: NCEP/NCAR再分析资料
Temporal Coverage: 从1948年1月1日起，一天4次/逐日/逐月数据
-----> 我们下载的是monthly也就是逐月数据
Spatial Coverage: 分辨率为2.5×2.5°的格点，0.0E~357.5E, 90.0N~90.0S
Levels: 垂直分为17层，分别是1000hPa, 925hPa......以此类推
那除了这些信息，更多的数据信息是存储在.nc文件本身里的，对于这些数据，我们可以直接通过ncl进行查看。
数据信息的查看
1、ncl_filedump命令(未进入ncl)：ncl_filedump文档
打开终端，输入ncl_filedump hgt.mon.mean.nc，注意数据文件的路径！！要么先cd到当前目录下（如下图），要么在.nc文件名前加入相对路径！
Fig.3.png
2、ncl导入数据后再查看：
打开终端，输入ncl按回车进入ncl命令行模式中。输入：
f = addfile("..../hgt.mon.mean.nc", "r") ; 将数据读取并写入f这个变量中
printVarSummary(f) ; 输出f变量的主要信息
print(f) ; 输出f变量的所有信息
可以对比一下两种查看方式的异同点~
Fig.4.对比.png
数据信息含义
由以上两种方式输出的数据信息如下（以注释形式在代码中标出，每行;之后的内容为注释）：
Variable: f
Type: file
filename: hgt.mon.mean
path: hgt.mon.mean.nc
file global attributes:
description : Data from NCEP initialized reanalysis (4x/day). These are interpolated to pressure surfaces from model (sigma) surfaces.
platform : Model
Conventions : COARDS
NCO : 20121012
history : Created by NOAA-CIRES Climate Diagnostics Center (SAC) from the NCEP
reanalysis data set on 07/07/97 by calc.mon.mean.year.f using
/Datasets/nmc.reanalysis.derived/pressure/hgt.mon.mean.nc
from /Datasets/nmc.reanalysis/pressure/hgt.79.nc to hgt.95.nc
Converted to chunked, deflated non-packed NetCDF4 2014/09
title : monthly mean hgt from the NCEP Reanalysis
References : http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.html
dataset_title : NCEP-NCAR Reanalysis 1
dimensions:
level = 17 ; level这个维度有17层
lat = 73 ; 纬度73个
lon = 144 ; 精度144个
time = 832 // unlimited
variables:
float level ( level ) ; 变量名:level
units : millibar ; 单位:hPa
long_name : Level
positive : down
GRIB_id : 100
GRIB_name : hPa
actual_range : ( 1000, 10 )
axis : Z
float lat ( lat ) ; 变量名:lat
units : degrees_north ; 单位: 北纬
actual_range : ( 90, -90 ) ; 值的范围
long_name : Latitude
standard_name : latitude
axis : Y
float lon ( lon )
units : degrees_east
long_name : Longitude
actual_range : ( 0, 357.5 )
standard_name : longitude
axis : X
double time ( time )
long_name : Time
delta_t : 0000-01-00 00:00:00
avg_period : 0000-01-00 00:00:00
prev_avg_period : 0000-00-01 00:00:00
standard_name : time
axis : T
units : hours since 1800-01-01 00:00:0.0
actual_range : ( 1297320, 1904352 )
float hgt ( time, level, lat, lon ) ; 变量名:hgt
long_name : Monthly mean geopotential height
valid_range : ( -700, 35000 )
units : m
add_offset : 0
scale_factor : 1
missing_value : -9.96921e+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 : NCEP Reanalysis Derived Products
actual_range : ( -354.4583, 32321.1 )
_FillValue : -9.96921e+36
可以看到，其中hgt这个变量就是我们要用的，在某等压面下的位势数据啦。
float hgt ( time, level, lat, lon )说明hgt这个变量是个四维(832*17*73*144)的变量。
首先将这个变量从f中抽取出来，即输入：
ncl 0> f = addfile("hgt.mon.mean.nc", "r")
ncl 1> hgt = f->hgt
ncl 2>
如果不放心还可以使用printVarSummary()函数输出一下变量
ncl 0> f = addfile("hgt.mon.mean.nc", "r")
ncl 1> hgt = f->hgt
ncl 2> printVarSummary(hgt)
Variable: hgt
Type: float
Total Size: 594726912 bytes
148681728 values
Number of Dimensions: 4
Dimensions and sizes: [time | 832] x [level | 17] x [lat | 73] x [lon | 144]
Coordinates:
time: [1297320..1904352]
level: [1000..10]
lat: [90..-90]
lon: [ 0..357.5]
Number Of Attributes: 17
long_name : Monthly mean geopotential height
valid_range : ( -700, 35000 )
units : m
add_offset : 0
scale_factor : 1
missing_value : -9.96921e+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 : NCEP Reanalysis Derived Products
actual_range : ( -354.4583, 32321.1 )
_FillValue : -9.96921e+36
ncl 3>
数据的切片和整合
命题为：【画出1948年以来500hPa高度上的平均等位势线】
目前hgt这个变量time有832维，level也是包括17个等压面的数据，我们需要：
1、将时间维度做平均
2、切片只留500hPa这个level的数据
1、时间维度平均
ncl 3> mean_hgt = dim_avg_n_Wrap(hgt, 0) ; 0表示对hgt的第一维做平均
输出的mean_hgt如下：
ncl 4> printVarSummary(mean_hgt)
Variable: mean_hgt
Type: float
Total Size: 714816 bytes
178704 values
Number of Dimensions: 3
Dimensions and sizes: [level | 17] x [lat | 73] x [lon | 144]
Coordinates:
level: [1000..10]
lat: [90..-90]
lon: [ 0..357.5]
Number Of Attributes: 18
_FillValue : -9.96921e+36
long_name : Monthly mean geopotential height
valid_range : ( -700, 35000 )
units : m
add_offset : 0
scale_factor : 1
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 : NCEP Reanalysis Derived Products
actual_range : ( -354.4583, 32321.1 )
missing_value : -9.96921e+36
average_op_ncl : dim_avg_n over dimension(s): time
2、切片只留level为500hPa的数据
切片操作我想大家都会，那么问题来了，500hPa是第几层？
level是从高到底存的，还是从低到高？
float hgt ( time, level, lat, lon )这行代码同时也意味着，level这个维度，对应的坐标是存在level这个变量中的，其顺序和level变量中的存储顺序是相同的。也就是说，我们看level这个变量里是什么样的顺序就好了。
ncl 5> print(f->level)
Variable: level (file variable)
Type: float
Total Size: 68 bytes
17 values
Number of Dimensions: 1
Dimensions and sizes: [level | 17]
Coordinates:
level: [1000..10]
Number Of Attributes: 7
units : millibar
long_name : Level
positive : down
GRIB_id : 100
GRIB_name : hPa
actual_range : ( 1000, 10 )
axis : Z
(0) 1000
(1) 925
(2) 850
(3) 700
(4) 600
(5) 500 ; 看这里看这里看这里！！！
(6) 400
(7) 300
(8) 250
(9) 200
(10) 150
(11) 100
(12) 70
(13) 50
(14) 30
(15) 20
(16) 10
看到了吧，500hPa是第6个值，但是由于序号从0开始标，所以index=5.
ncl 6> mean_500_hgt = mean_hgt(5, :, :)
ncl 7> print(mean_500_hgt)
Variable: mean_500_hgt
Type: float
Total Size: 42048 bytes
10512 values
Number of Dimensions: 2
Dimensions and sizes: [lat | 73] x [lon | 144]
Coordinates:
lat: [90..-90]
lon: [ 0..357.5]
Number Of Attributes: 19
level : 500
average_op_ncl : dim_avg_n over dimension(s): time
missing_value : -9.96921e+36
actual_range : ( -354.4583, 32321.1 )
dataset : NCEP Reanalysis Derived Products
parent_stat : Other
statistic : Mean
level_desc : Multiple levels
var_desc : Geopotential height
GRIB_name : HGT
GRIB_id : 7
least_significant_digit : 0
precision : 0
scale_factor : 1
add_offset : 0
units : m
valid_range : ( -700, 35000 )
long_name : Monthly mean geopotential height
_FillValue : -9.96921e+36
现在这个变量是一个二维变量了~ 我们已经可以将其画在地图上
画图设置
1、设置工作台
wks = gsn_open_wks("png", "500hPa mean hgt")
这是最简单的工作台设置方式，
第一个参数表示生成的格式，"png"表示png格式的图片
第二个参数是生成的文件名称
这样设置出来以后意味着会生成一个500hPa mean hgt.png的文件，一会儿画出来图就去相应的路径下找这个文件就好啦~
2、设置图片属性
这一部分对于强迫症患者来说真的是最最最恶心的部分了，NCL语言设计中最最最恶心的地方就是在这里。图片属性中的每一个值都有对应的参数名称，你想修改任意一个点都必须找到相应的参数，然后给予设置，不然就都是默认值！！
以下是最简单的图片属性设置:
res = True ; resource变量为True
res@cnFillOn = True ; 是否填充颜色
res@cnLinesOn = False ; 是否绘制等值线
res@cnLineLabelsOn = False ; 是否显示等值线上的值标签
; cn这一溜表示这些参数是contour图的属性
画图
设置好了工作台和图片属性，画图其实就是一行代码的事。
plot = gsn_csm_contour_map(wks, mean_500_hgt, res)
gsn_csm_contour_map是一个绘图函数，_map表示绘制在地图上，_contour表示画图的种类，除了contour还有xy图等等等好多，以后会讲。
第一个参数，填入工作台的变量。
第二个参数，是你要绘制的变量。
第三个参数，是存储图片属性的变量。
绘制出来的图片：
Fig.5.500hPa平均位势高度.png
....看吧，位势高度场的年平均值……画出来真的是毫无意义啊……有兴趣的同学可以画逐月的值，或者是夏季平均值~
结束语
只是为了让大家快速地对NCL整体绘图流程有一个模糊的概念，我相信中间有很多地方大家有疑惑的，我会努力尽快整理出来~
大家对文章有哪些建议或者意见，还请直接留言评论~
有哪些迫切想要了解的点，也可以告诉我，我可以提前梳理出来。
附：完整代码
; ************************************
; Chapter.2.开始一个简单的NCL画图程序
; ************************************
; 直接打开一个编辑器，另存为.ncl的文件就行了
; 自己回顾一下完整的.ncl代码要怎么跑！
; *************************************
begin
f = addfile("hgt.mon.mean.nc","r")
; printVarSummary(f)
hgt = f->hgt
mean_hgt = dim_avg_n_Wrap(hgt, 0)
; printVarSummary(mean_hgt)
;print(f->level)
mean_500_hgt = mean_hgt(5, :, :)
; ------ set the work station ------
wks = gsn_open_wks("png", "500hPa mean hgt")
; ------ set the resources ------
res = True
res@cnFillOn = True
res@cnLinesOn = False
res@cnLineLabelsOn = False
plot= gsn_csm_contour_map(wks, mean_500_hgt, res)
end
网友评论