首页 > 代码库 > NCL 学习笔记
NCL 学习笔记
这里以绘制气温分布图为例，效果如下图：
这里几点说明：
1.ncl不支持中文显示，所有文字都是英文，但是支持很多样式的字体，参考网址http://www.ncl.ucar.edu/Document/Graphics/font_tables.shtml
2.图下方的labelbar只能在图的周围，不能放置在图内。要想显示图下方的图例，就要使用legend而不是labelbar了。
使用NCL脚本绘制一张如上图所示的png图片主要分为以下几个步骤
一、读取各站点的气温数据。
二、将站点数据使用各种差值函数转换成格点数据。
三、使用源对地图进行基本设置
四、使用源对等值线填充进行基本设置
五、使用源对labelbar进行基本设置
六、生成png图片
接下来将按照这几个步骤，详细介绍。
一、读取各站点的气温数据
NCL支持的数据格式主要有netCDF文件（.nc .cdf）、HDF4（.hd .hdf）、HDF4-EOS(.hdfeos)、GRID-1/GRIB-2(.grb.grib)、CCMHistory Tape(.ecm),除此之外呢，它支持二进制文件和ascii文件，这两者是我们最熟悉的。这里我们使用ascii文件，更多文件读取方式参考http://www.ncl.ucar.edu/Applications/list_io.shtml
为了批量生成产品图片，需要配置文件设置数据来源以及图片生成后存放位置。config.txt文件如下：
One Hour of Temperature 2010111502
./t1/
/root/WorkSpace/MICAPS_surface/t1/
10111502.000
第一行是标题
第二行是输出png图路径
第三行是输入数据文件路径
第四行是数据文件名
在NCL脚本（temperature.ncl）中使用以下几行代码就可以了
filepath = "./config.txt" ;参数文件路径
argu = asciiread(filepath,-1,"string") ;以字符串形式读取参数文件入数组argu
lines = asciiread(argu(2)+argu(3),-1,"string") ;以字符串形式读取数据文件入数组lines
station = stringtofloat(str_get_field(lines(3::),1," ")) ;从数组lines中获取站号
lon = stringtofloat(str_get_field(lines(3::),2," ")) ;从数组lines中获取经度值lon
lat = stringtofloat(str_get_field(lines(3::),3," ")) ;从数组lines中获取纬度值lat
height = stringtofloat(str_get_field(lines(3::),4," ")) ;从数组lines中获取海拔高度
R = stringtofloat(str_get_field(lines(3::),5," ")) ;从数组lines中获取站点数据值
由于数据文件10111502.000的前3行是文件头，不包含数据，因此lines从第三行开始读取数据。注意NCL中注释使用“;”而且只能注释单行。
这样所有数据就保存到各个变量中啦！每个变量都是一个一维数组。
二、将站点数据使用各种差值函数转换成格点数据。
接下来使用插值函数，NCL中提供了很多差值函数，如Cressman插值（obj_anal_ic_deprecated）、三次样条差值（caa1）、反距离权重差值（dsgrid2）、最邻近点产值（natgrid）等等，其他插值函数参考网址http://www.ncl.ucar.edu/Document/Functions/interp.shtml
注意使用这些函数的时候要在文件头部包含contributed.ncl，即：
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
各个差值函数的接口不同，需要提前定义的变量也会有所差异，这里使用Cressman插值。
首先创建存放差值后生成数据的数组 根据亚洲截图的经纬度而定北纬17-57，东经72-138这个矩形框内插值
olon = new(66,"float");
olat = new(40,"float");
data1 = new((/40,66/),"float")
do i=0,65
olon(i) =72+i
end do
do l=0,39
olat(l) = 17+l
end do
接下来设置数组属性，为了符合netcdf规定的数据格式，使函数能够识别经纬度
olon!0 = "lon"
olon@long_name = "lon"
olon@units = "degrees-east"
olon&lon = olon
olat!0 = "lat"
olat@long_name = "lat"
olat@units = "degrees_north"
olat&lat = olat
最后调用插值函数
R@_FillValue = 999999.000000
rscan = (/10,5,3/) ;连续的有效半径大小，最大为10，依次递减
data1 = obj_anal_ic_deprecated(lon,lat,R,olon,olat,rscan, False) ;Creanm插值
三、使用源对地图进行基本设置
首先创建一个自定义的colormap配色方案，并且创建一个工作空间，
cmap = (/ \
(/ 255./255, 255./255, 255./255 /), \ ; 0 - White background.
(/ 0./255 , 0./255 , 0./255 /), \ ; 1 - Black foreground.
(/ 255./255, 0./255 , 0./255 /), \ ; 2 - Red.
(/ 0./255 , 0./255 , 255./255 /), \ ; 3 - Blue.
(/ 164./255, 244./255, 131./255 /), \ ; 4 - Ocean Blue.
(/ 0./255 , 0./255 , 255./255 /), \ ; 5 - Bar 1
(/ 0./255 , 153./255, 255./255 /), \ ; 6 - Bar 2
(/ 0./255, 153./255, 153./255 /), \ ; 7 - Bar 3
(/ 0./255 , 255./255, 0./255 /), \ ; 8 - Bar 4
(/ 255./255, 255./255 , 102./255 /), \ ; 9 - Bar 5
(/ 255./255, 153./255 , 102./255 /), \ ; 10 - Bar 6
(/ 255./255, 0./255 , 255./255 /) \ ; 11 - Bar 7
/)
wks_type = "png"
wks = gsn_open_wks(wks_type,argu(1)+argu(3)) ; Open a workstation and.
gsn_define_colormap(wks,cmap) ; define a different colormap.
这里创建了png的工作空间，NCL还支持X11、PS、NCGM、PDF、NEWPDF等。
接下来设置地图属性
res = True
res@gsnAddCyclic = False ;由于我们的数据不是循环地球一周的，因此必须把这个置否
res@mpDataSetName = "Earth..4" ; This new database contains
; divisions for other countries.
res@mpDataBaseVersion = "MediumRes" ; Medium resolution database
res@mpOutlineOn = True ; Turn on map outlines
res@mpOutlineSpecifiers = (/"China:states","Taiwan"/) ;China:states
中国地图包含在Earth..4这个地图库中，将边界区域设置为中国行政区域和台湾，在台湾问题这一点上比较郁闷，中国地图里没有台湾，激起了我这个爱国主义青年的强烈愤慨。
地图选好了，该把区域缩小到中国范围内了，这里和上面的插值范围有些出入，只是显示需要，没有实质联系。
res@mpMinLatF = 17 ; Asia limits
res@mpMaxLatF = 55
res@mpMinLonF = 72
res@mpMaxLonF = 136
你还可以使用这两行代码来加粗边界线。
res@mpGeophysicalLineThicknessF= 2. ; double the thickness of geophysical boundaries
res@mpNationalLineThicknessF= 2. ; double the thickness of national boundaries
默认的底图投影方式是等经纬度投影，画出来的中国地图比较扁，我们常看到的中国地图，投影方式是兰伯特投影，因此需要对投影方式进行修改
res@mpProjection = "LambertConformal" ;兰伯特投影
res@mpLambertMeridianF = 110.0
res@mpLimitMode = "LatLon"
res@mpLambertParallel1F = .001 ;Default: .001 ;可以自己改一改，看看投影有什么不同，挺有趣的
res@mpLambertParallel2F = 89.999 ;Default: 89.999
最后将填充区域设定在中国行政区域图之内，如果使用默认效果，等值线会对整个矩形区域填充颜色，因此需要去掉中国边境范围外的填充颜色
res@mpAreaMaskingOn = True ;使能填充覆盖
res@mpMaskAreaSpecifiers = (/"China:states","Taiwan"/) ;China:states
res@mpOceanFillColor = 0 ;用白色填充海洋 0是colormap的索引值
res@mpInlandWaterFillColor = 0 ;用白色填充内陆湖水
四、使用源对等值线填充进行基本设置
首先使能等值线填充，不显示各填充颜色之间的黑色等值线，不显示等值线上标示的数值，在绘制其他要素前先绘制等值线，对源的设置如下：
res@cnFillOn = True
res@cnLinesOn = False ;等值线不显示
res@cnLineLabelsOn = False
res@cnFillDrawOrder = "PreDraw" ; draw contours first
使用特定的Level值和配色方案对等值区间进行设置
res@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels
res@cnLevels = (/-30.,-20.,-10.,0.,10.,20./) ; set levels
res@cnFillColors = (/5,6,7,8,9,10,11/) ; set the colors to be used
五、使用源对labelbar进行基本设置
res@lbLabelBarOn = True ;LabelBar显示
res@lbLabelStrings = (/"-30:S:o:N","-20:S:o:N","-10:S:o:N","0:S:o:N","10:S:o:N","20:S:o:N"/)
其中“:S:o:N”表示摄氏度符号
你也可以使用以下这行代码，将labelbar放置在图片右侧
res@lbOrientation = "vertical" ; vertical label bar
六、生成png图片
在最后生成图片前，先要对整个view进行调整
res@vpXF = 0 ;左边距
res@vpYF = 0.95 ;上边距
res@vpWidthF = 1.0 ; height and width of plot
res@vpHeightF = 0.8
然后再设置图片标题
res@tiMainFont = "helvetica"
res@tiMainOffsetYF = 0.02 ;set place for main title along Y,offset
res@tiMainFontHeightF = 0.02 ;set main title font size
res@tiMainString = argu(0)
最后显示
map = gsn_csm_contour_map(wks,data1,res)
别忘了在整个代码前后增加begin（代码）end这种固定格式。
这样这张图片就顺利生成啦！
七、高级进阶
注意到样例图片上各个省都是有汉语拼音名称的，NCL可不会自动帮你把它填上去，它的实现需要我们自己写代码，具体实现方法也不难，只需要提供各个省的经纬度，以及需要填充的字符串即可。
因为我们是在原有的map基础上再次绘图，所以需要现对map的res源进行设置，不让它绘制后就立马显示，增加代码res@gsnFrame = False ; don‘t advance frame yets
这样我们就可以在添加完各省名称后统一显示。
另外创建一个tres源用于显示文本，
tres = True ; text mods desired
tres@txFontHeightF = 0.01 ; make smaller
plat = (/31.69,39.76,29.37,26.01,36.46,\
22.96,23.69,26.4,19.18,38.19,47.04,\
34.64,22.03,30.99,28.22,33.53,27.67,43.82,\
41.54,40.74,38.0,35.14,34.77,\
36.5,31.06,37.49,30.83,38.95,\
40.11,30.75,24.26,29.3,24.09/)
plon = (/117.2,116.3,106.5,118.1,101.7,\
113.2,109.1,106.7,109.8,115.2,129,\
113.7,113.5,112.7,111.3,119.2,115.5,125.9,\
123.1,109.1,107.2,96.53,109.3,\
117.0,121.4,112.2,102.3,117.2,\
87.65,90.99,101.2,119.2,120.9/)
do i=0,32
gsn_text(wks,map,province(i),plon(i),plat(i),tres)
end do
最后调用
frame(wks)
显示叠加后的图片。
NCL 学习笔记