这里以绘制气温分布图为例,效果如下图:
这里几点说明:
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)、CCMHistoryTape(.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 - Bar4
(/ 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 - Bar7
/)
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 newdatabase 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 contoursfirst
使用特定的Level值和配色方案对等值区间进行设置
res@cnLevelSelectionMode ="ExplicitLevels"; setexplicit 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 beused
五、使用源对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 alongY,offset
res@tiMainFontHeightF= 0.02 ;set main title fontsize
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)
显示叠加后的图片。
终于大功告成啦!