NCL学习笔记(实战篇)

这里以绘制气温分布图为例，效果如下图：
X 这里几点说明：
1.ncl不支持中文显示，所有文字都是英文，但是支持很多样式的字体，参考
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文件，更多文件读取方式参考/Applications/list_io.shtml
为了批量生成产品图片，需要配置文件设置数据来源以及图片生成后存放位置。

config.txt文件如下：
One Hour of Temperature2010111502
./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）等等，其他插值函数参考网址/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 ;由于我们的数据不是循环地球一周的，因此必须把这个置否
中国地图包含在Earth..4这个地图库中，将边界区域设置为中国行政区域和台湾，在台湾问题这一点上比较郁闷，中国地图里没有台湾，激起了我这个爱国主义青年的强烈愤慨。

地图选好了，该把区域缩小到中国范围内了，这里和上面的插值范围有些出入，只是显示需要，没有实质联系。

你还可以使用这两行代码来加粗边界线。