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)、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
NCL学习笔记(实战篇)
(/ 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)

显示叠加后的图片。

终于大功告成啦!

  

爱华网本文地址 » http://www.aihuau.com/a/25101018/368050.html

更多阅读

中小学教师职业道德规范学习笔记

中小学教师职业道德规范学习笔记一、依法执教。学习和宣传马列主义、毛泽东思想和邓--同志建设有中国特色社会主义理论,拥护党的基本路线,全面贯彻国家教育方针,自觉遵守《教师法》等法律法规,在教育教学中同党和国家的方针政策保持一

伤寒学习笔记(下

伤寒学习笔记---14,桂枝用量问题?学习仲景药法,涉及用量问题。深入考证,固然不必,简单了解,实属必要。仅以桂枝为例,整理如下:基本资料:《现代中医药应用与研究大系》1985年版,简称《大系》。药法:凡例里说:“将其折合今之用量---,并结合笔者临床

关于闪光灯TTL的学习笔记

关于闪光灯TTL的学习笔记(转)俺注:一直对加闪光灯的拍摄比较疑惑,要知道并不是加开个闪光就一了百了的。特别是现在的专业闪灯,有好几种模式,都需要配合相机的测光使用。今天无意之间看到一篇解疑的文章,赶紧收集下来仔细学习。想学习下闪

声明:《NCL学习笔记实战篇》为网友寄与情人分享!如侵犯到您的合法权益请联系我们删除