三菱欧蓝德火花塞:[转载]NCL学习笔记（实战篇）

[转载]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)、CCMHistory Tape(.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 - 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      ;由于我们的数据不是循环地球一周的，因此必须把这个置否

  res@mpDataSetName         = "Earth..4"   ; This new database 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 contours first   

    使用特定的Level值和配色方案对等值区间进行设置

  res@cnLevelSelectionMode = "ExplicitLevels"       ; set explicit 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 be used 

五、使用源对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 along Y,offset
  res@tiMainFontHeightF     = 0.02   ;set main title font size
  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)

显示叠加后的图片。

终于大功告成啦！