【zz】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
           (/ 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)

显示叠加后的图片。