這里以繪制氣溫分布圖為例,效果如下圖:
這里幾點(diǎn)說明:
1.ncl不支持中文顯示,所有文字都是英文,但是支持很多樣式的字體,參考網(wǎng)址http://www.ncl.ucar.edu/Document/Graphics/font_tables.shtml
2.圖下方的labelbar只能在圖的周圍,不能放置在圖內(nèi)。要想顯示圖下方的圖例,就要使用legend而不是labelbar了。
使用NCL腳本繪制一張如上圖所示的png圖片主要分為以下幾個(gè)步驟
一、讀取各站點(diǎn)的氣溫?cái)?shù)據(jù)。
二、將站點(diǎn)數(shù)據(jù)使用各種差值函數(shù)轉(zhuǎn)換成格點(diǎn)數(shù)據(jù)。
三、使用源對地圖進(jìn)行基本設(shè)置
四、使用源對等值線填充進(jìn)行基本設(shè)置
五、使用源對labelbar進(jìn)行基本設(shè)置
六、生成png圖片
接下來將按照這幾個(gè)步驟,詳細(xì)介紹。
一、讀取各站點(diǎn)的氣溫?cái)?shù)據(jù)
NCL支持的數(shù)據(jù)格式主要有netCDF文件(.nc.cdf)、HDF4(.hd.hdf)、HDF4-EOS(.hdfeos)、GRID-1/GRIB-2(.grb.grib)、CCMHistoryTape(.ecm),除此之外呢,它支持二進(jìn)制文件和ascii文件,這兩者是我們最熟悉的。這里我們使用ascii文件,更多文件讀取方式參考http://www.ncl.ucar.edu/Applications/list_io.shtml
為了批量生成產(chǎn)品圖片,需要配置文件設(shè)置數(shù)據(jù)來源以及圖片生成后存放位置。config.txt文件如下:
One Hour of Temperature 2010111502
./t1/
/root/WorkSpace/MICAPS_surface/t1/
10111502.000
第一行是標(biāo)題
第二行是輸出png圖路徑
第三行是輸入數(shù)據(jù)文件路徑
第四行是數(shù)據(jù)文件名
在NCL腳本(temperature.ncl)中使用以下幾行代碼就可以了
filepath ="./config.txt" ;參數(shù)文件路徑
argu = asciiread(filepath,-1,"string");以字符串形式讀取參數(shù)文件入數(shù)組argu
lines = asciiread(argu(2)+argu(3),-1,"string");以字符串形式讀取數(shù)據(jù)文件入數(shù)組lines
station =stringtofloat(str_get_field(lines(3::),1," ")) ;從數(shù)組lines中獲取站號
lon =stringtofloat(str_get_field(lines(3::),2," "));從數(shù)組lines中獲取經(jīng)度值lon
lat =stringtofloat(str_get_field(lines(3::),3," "));從數(shù)組lines中獲取緯度值lat
height =stringtofloat(str_get_field(lines(3::),4," ")) ;從數(shù)組lines中獲取海拔高度
R =stringtofloat(str_get_field(lines(3::),5," "));從數(shù)組lines中獲取站點(diǎn)數(shù)據(jù)值
由于數(shù)據(jù)文件10111502.000的前3行是文件頭,不包含數(shù)據(jù),因此lines從第三行開始讀取數(shù)據(jù)。注意NCL中注釋使用“;”而且只能注釋單行。
這樣所有數(shù)據(jù)就保存到各個(gè)變量中啦!每個(gè)變量都是一個(gè)一維數(shù)組。
二、將站點(diǎn)數(shù)據(jù)使用各種差值函數(shù)轉(zhuǎn)換成格點(diǎn)數(shù)據(jù)。
接下來使用插值函數(shù),NCL中提供了很多差值函數(shù),如Cressman插值(obj_anal_ic_deprecated)、三次樣條差值(caa1)、反距離權(quán)重差值(dsgrid2)、最鄰近點(diǎn)產(chǎn)值(natgrid)等等,其他插值函數(shù)參考網(wǎng)址http://www.ncl.ucar.edu/Document/Functions/interp.shtml
注意使用這些函數(shù)的時(shí)候要在文件頭部包含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"
各個(gè)差值函數(shù)的接口不同,需要提前定義的變量也會(huì)有所差異,這里使用Cressman插值。
首先創(chuàng)建存放差值后生成數(shù)據(jù)的數(shù)組 根據(jù)亞洲截圖的經(jīng)緯度而定北緯17-57,東經(jīng)72-138這個(gè)矩形框內(nèi)插值
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
接下來設(shè)置數(shù)組屬性,為了符合netcdf規(guī)定的數(shù)據(jù)格式,使函數(shù)能夠識別經(jīng)緯度
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
最后調(diào)用插值函數(shù)
R@_FillValue =999999.000000
rscan =(/10,5,3/);連續(xù)的有效半徑大小,最大為10,依次遞減
data1 =obj_anal_ic_deprecated(lon,lat,R,olon,olat,rscan,False) ;Creanm插值
三、使用源對地圖進(jìn)行基本設(shè)置
首先創(chuàng)建一個(gè)自定義的colormap配色方案,并且創(chuàng)建一個(gè)工作空間,
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.
這里創(chuàng)建了png的工作空間,NCL還支持X11、PS、NCGM、PDF、NEWPDF等。
接下來設(shè)置地圖屬性
res = True
res@gsnAddCyclic =False;由于我們的數(shù)據(jù)不是循環(huán)地球一周的,因此必須把這個(gè)置否
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這個(gè)地圖庫中,將邊界區(qū)域設(shè)置為中國行政區(qū)域和臺灣,在臺灣問題這一點(diǎn)上比較郁悶,中國地圖里沒有臺灣,激起了我這個(gè)愛國主義青年的強(qiáng)烈憤慨。
地圖選好了,該把區(qū)域縮小到中國范圍內(nèi)了,這里和上面的插值范圍有些出入,只是顯示需要,沒有實(shí)質(zhì)聯(lián)系。
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
默認(rèn)的底圖投影方式是等經(jīng)緯度投影,畫出來的中國地圖比較扁,我們??吹降闹袊貓D,投影方式是蘭伯特投影,因此需要對投影方式進(jìn)行修改
res@mpProjection ="LambertConformal" ;蘭伯特投影
res@mpLambertMeridianF =110.0
res@mpLimitMode = "LatLon"
res@mpLambertParallel1F =.001;Default:.001;可以自己改一改,看看投影有什么不同,挺有趣的
res@mpLambertParallel2F =89.999;Default: 89.999
最后將填充區(qū)域設(shè)定在中國行政區(qū)域圖之內(nèi),如果使用默認(rèn)效果,等值線會(huì)對整個(gè)矩形區(qū)域填充顏色,因此需要去掉中國邊境范圍外的填充顏色
res@mpAreaMaskingOn =True ;使能填充覆蓋
res@mpMaskAreaSpecifiers =(/"China:states","Taiwan"/);China:states
res@mpOceanFillColor =0;用白色填充海洋 0是colormap的索引值
res@mpInlandWaterFillColor= 0 ;用白色填充內(nèi)陸湖水
四、使用源對等值線填充進(jìn)行基本設(shè)置
首先使能等值線填充,不顯示各填充顏色之間的黑色等值線,不顯示等值線上標(biāo)示的數(shù)值,在繪制其他要素前先繪制等值線,對源的設(shè)置如下:
res@cnFillOn= True
res@cnLinesOn=False;等值線不顯示
res@cnLineLabelsOn = False
res@cnFillDrawOrder="PreDraw"; draw contoursfirst
使用特定的Level值和配色方案對等值區(qū)間進(jìn)行設(shè)置
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進(jìn)行基本設(shè)置
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放置在圖片右側(cè)
res@lbOrientation="vertical"; vertical label bar
六、生成png圖片
在最后生成圖片前,先要對整個(gè)view進(jìn)行調(diào)整
res@vpXF =0;左邊距
res@vpYF =0.95;上邊距
res@vpWidthF =1.0; height and width of plot
res@vpHeightF= 0.8
然后再設(shè)置圖片標(biāo)題
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)
別忘了在整個(gè)代碼前后增加begin(代碼)end這種固定格式。
這樣這張圖片就順利生成啦!
七、高級進(jìn)階
注意到樣例圖片上各個(gè)省都是有漢語拼音名稱的,NCL可不會(huì)自動(dòng)幫你把它填上去,它的實(shí)現(xiàn)需要我們自己寫代碼,具體實(shí)現(xiàn)方法也不難,只需要提供各個(gè)省的經(jīng)緯度,以及需要填充的字符串即可。
因?yàn)槲覀兪窃谠械膍ap基礎(chǔ)上再次繪圖,所以需要現(xiàn)對map的res源進(jìn)行設(shè)置,不讓它繪制后就立馬顯示,增加代碼res@gsnFrame=False; don't advance frame yets
這樣我們就可以在添加完各省名稱后統(tǒng)一顯示。
另外創(chuàng)建一個(gè)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
最后調(diào)用
frame(wks)
顯示疊加后的圖片。
終于大功告成啦!
愛華網(wǎng)



