為了更好地對(duì)NCL有一個(gè)整體的認(rèn)知葵硕,首先竞穷,舉一個(gè)最簡(jiǎn)單的例子梳庆。
畫出1948年以來(lái)500hPa高度上的平均等位勢(shì)線暖途。
數(shù)據(jù)可從以下網(wǎng)址中獲取:https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html
有了數(shù)據(jù)膏执,知道了要畫什么驻售,首先你要做的是:知道數(shù)據(jù)的存儲(chǔ)格式,也就是知道數(shù)據(jù)長(zhǎng)什么樣子更米。
像我們剛剛下載的數(shù)據(jù)芋浮,屬于比較正規(guī)的數(shù)據(jù),其網(wǎng)站上面有一些說(shuō)明:
Brief Description: NCEP/NCAR再分析資料
Temporal Coverage: 從1948年1月1日起壳快,一天4次/逐日/逐月數(shù)據(jù)
-----> 我們下載的是monthly也就是逐月數(shù)據(jù)
Spatial Coverage: 分辨率為2.5×2.5°的格點(diǎn)纸巷,0.0E~357.5E, 90.0N~90.0S
Levels: 垂直分為17層,分別是1000hPa, 925hPa......以此類推
那除了這些信息眶痰,更多的數(shù)據(jù)信息是存儲(chǔ)在.nc文件本身里的瘤旨,對(duì)于這些數(shù)據(jù),我們可以直接通過ncl進(jìn)行查看竖伯。
數(shù)據(jù)信息的查看
1存哲、ncl_filedump
命令(未進(jìn)入ncl):ncl_filedump文檔
打開終端,輸入ncl_filedump hgt.mon.mean.nc
七婴,注意數(shù)據(jù)文件的路徑K钔怠!要么先cd
到當(dāng)前目錄下(如下圖)打厘,要么在.nc文件名前加入相對(duì)路徑修肠!
2、ncl導(dǎo)入數(shù)據(jù)后再查看:
打開終端户盯,輸入ncl
按回車進(jìn)入ncl命令行模式中嵌施。輸入:
f = addfile("..../hgt.mon.mean.nc", "r") ; 將數(shù)據(jù)讀取并寫入f這個(gè)變量中
printVarSummary(f) ; 輸出f變量的主要信息
print(f) ; 輸出f變量的所有信息
可以對(duì)比一下兩種查看方式的異同點(diǎn)~
數(shù)據(jù)信息含義
由以上兩種方式輸出的數(shù)據(jù)信息如下(以注釋形式在代碼中標(biāo)出,每行;之后的內(nèi)容為注釋):
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這個(gè)維度有17層
lat = 73 ; 緯度73個(gè)
lon = 144 ; 精度144個(gè)
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 ; 缺失值填補(bǔ)為這個(gè)值
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
這個(gè)變量就是我們要用的吗伤,在某等壓面下的位勢(shì)數(shù)據(jù)啦。
float hgt ( time, level, lat, lon )
說(shuō)明hgt這個(gè)變量是個(gè)四維(832*17*73*144)的變量硫眨。
首先將這個(gè)變量從f中抽取出來(lái)足淆,即輸入:
ncl 0> f = addfile("hgt.mon.mean.nc", "r")
ncl 1> hgt = f->hgt
ncl 2>
如果不放心還可以使用printVarSummary()
函數(shù)輸出一下變量
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>
數(shù)據(jù)的切片和整合
命題為:【畫出1948年以來(lái)500hPa高度上的平均等位勢(shì)線】
目前hgt這個(gè)變量time有832維,level也是包括17個(gè)等壓面的數(shù)據(jù)礁阁,我們需要:
1巧号、將時(shí)間維度做平均
2、切片只留500hPa這個(gè)level的數(shù)據(jù)
1氮兵、時(shí)間維度平均
ncl 3> mean_hgt = dim_avg_n_Wrap(hgt, 0) ; 0表示對(duì)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的數(shù)據(jù)
切片操作我想大家都會(huì),那么問題來(lái)了泣栈,500hPa是第幾層卜高?
level是從高到底存的,還是從低到高南片?
float hgt ( time, level, lat, lon )
這行代碼同時(shí)也意味著掺涛,level這個(gè)維度,對(duì)應(yīng)的坐標(biāo)是存在level這個(gè)變量中的疼进,其順序和level變量中的存儲(chǔ)順序是相同的薪缆。也就是說(shuō),我們看level這個(gè)變量里是什么樣的順序就好了伞广。
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 ; 看這里看這里看這里<鹈薄L鄣纭!
(6) 400
(7) 300
(8) 250
(9) 200
(10) 150
(11) 100
(12) 70
(13) 50
(14) 30
(15) 20
(16) 10
看到了吧减拭,500hPa是第6個(gè)值蔽豺,但是由于序號(hào)從0開始標(biāo),所以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
現(xiàn)在這個(gè)變量是一個(gè)二維變量了~ 我們已經(jīng)可以將其畫在地圖上
畫圖設(shè)置
1拧粪、設(shè)置工作臺(tái)
wks = gsn_open_wks("png", "500hPa mean hgt")
這是最簡(jiǎn)單的工作臺(tái)設(shè)置方式修陡,
第一個(gè)參數(shù)表示生成的格式,"png"表示png格式的圖片
第二個(gè)參數(shù)是生成的文件名稱
這樣設(shè)置出來(lái)以后意味著會(huì)生成一個(gè)500hPa mean hgt.png
的文件可霎,一會(huì)兒畫出來(lái)圖就去相應(yīng)的路徑下找這個(gè)文件就好啦~
2魄鸦、設(shè)置圖片屬性
這一部分對(duì)于強(qiáng)迫癥患者來(lái)說(shuō)真的是最最最惡心的部分了,NCL語(yǔ)言設(shè)計(jì)中最最最惡心的地方就是在這里癣朗。圖片屬性中的每一個(gè)值都有對(duì)應(yīng)的參數(shù)名稱拾因,你想修改任意一個(gè)點(diǎn)都必須找到相應(yīng)的參數(shù),然后給予設(shè)置斯棒,不然就都是默認(rèn)值6苤隆!
以下是最簡(jiǎn)單的圖片屬性設(shè)置:
res = True ; resource變量為True
res@cnFillOn = True ; 是否填充顏色
res@cnLinesOn = False ; 是否繪制等值線
res@cnLineLabelsOn = False ; 是否顯示等值線上的值標(biāo)簽
; cn這一溜表示這些參數(shù)是contour圖的屬性
畫圖
設(shè)置好了工作臺(tái)和圖片屬性荣暮,畫圖其實(shí)就是一行代碼的事庭惜。
plot = gsn_csm_contour_map(wks, mean_500_hgt, res)
gsn_csm_contour_map
是一個(gè)繪圖函數(shù),_map表示繪制在地圖上穗酥,_contour表示畫圖的種類护赊,除了contour還有xy圖等等等好多,以后會(huì)講砾跃。
第一個(gè)參數(shù)骏啰,填入工作臺(tái)的變量。
第二個(gè)參數(shù)抽高,是你要繪制的變量判耕。
第三個(gè)參數(shù),是存儲(chǔ)圖片屬性的變量翘骂。
繪制出來(lái)的圖片:
....看吧壁熄,位勢(shì)高度場(chǎng)的年平均值……畫出來(lái)真的是毫無(wú)意義啊……有興趣的同學(xué)可以畫逐月的值,或者是夏季平均值~
結(jié)束語(yǔ)
只是為了讓大家快速地對(duì)NCL整體繪圖流程有一個(gè)模糊的概念碳竟,我相信中間有很多地方大家有疑惑的草丧,我會(huì)努力盡快整理出來(lái)~
大家對(duì)文章有哪些建議或者意見,還請(qǐng)直接留言評(píng)論~
有哪些迫切想要了解的點(diǎn)莹桅,也可以告訴我昌执,我可以提前梳理出來(lái)。
附:完整代碼
; ************************************
; Chapter.2.開始一個(gè)簡(jiǎn)單的NCL畫圖程序
; ************************************
; 直接打開一個(gè)編輯器,另存為.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