五月最后一周的NCL笔记

大家好！又是一个边骂NCL弱智一边写NCL的星期过去了！
我来汇报这星期新发现的bug！！！（。

【Resource类】

 如果你也想给你的Panel图加上一个Color Bar，那你或许大概有可能也会用到如下resource设置：

1
2
3
4
5
6
7
8
9
10
11
12
panel_res                            = True        

panel_res@gsnPanelDebug              = True    ;输出各子图的ndc坐标信息
panel_res@gsnPanelLabelBar           = True    ;打开panel的Label Bar
panel_res@lbOrientation              = "Horizontal"
panel_res@pmLabelBarOrthogonalPosF   = -0.015  ;把Label Bar往下挪一点
panel_res@pmLabelBarWidthF           = 0.4     ;设置Bar的宽度
panel_res@lbLabelFontHeightF         = 0.012   ;设置Bar的字体大小
panel_res@gsnPanelBottom             = 0.05    ;放不下了，给画布底下加长点儿

panel_res@lbBoxLinesOn               = True
panel_res@lbBoxLineColor             = "grey77" ;设置Bar中框线的颜色

中间比较神奇的就是虽然大部分Panel图资源都通过plot manager（pm）实现，但其方向和字体宽度等都仍然由label bar（lb）分类下的属性设置。另外不得不特别提一句的就是

  • gsnPanelDebug：一个强烈推荐panel图排版不如人意的人试着打开的开关，不细说，但是着实好用！
  • pmLabelBarOrthogonalPosF/pmLabelBarParallelPosF：前者设定垂直于label bar方向的ndc坐标偏移，后者设定平行于label bar方向的ncd坐标偏移，因而对于水平和垂直方向的label bar二者的方向是不同的。ndc坐标也是又一个NCL神奇设定，简而言之就是以左上角为(x=0.00,y=1.00)，然后y轴向下递减，x轴向右增加的这么一个坐标系。上述设置中涉及到0.xx的数字也大概是以这个坐标系为参考的…吧。
    （附：NDC坐标示意图）

【Function类】

dtrend, dtrend_n, dtrend_msg_n

 使用最小二乘法去除数据的线性趋势，其中最后一个允许数据中包含缺测值。

1
2
3
dtrend(y, return_info) - 对最右边一维去趋势
dtrend_n(y, return_info, dim) - 对dim指定的维去趋势
dtrend_msg_n(x, y, remove_mean, return_info, dim) - 同上一个，但允许缺测值

 其中return_info如果设置为True，返回值就会包括两个attribute，一个是去掉的趋势slope，一个是我也不知道用来干嘛的y_intercept，均为一维序列，如果要用它俩画图需要先reshape一下。最后一个函数里的第一个参数x是y需要去趋势那一维的坐标信息，一般直接取y&time(以时间维为例)即可。remove_mean正常情况下设为True。这系列函数都会保留坐标和atr信息，无需拷贝坐标信息。

mask

 “保护”一个变量中符合给定条件的部分，而把其他都设为缺测值。

 八百年搞不清楚的一个函数，每用一次都要查一次，常查常忘。此函数接收三个参数：

1
mask(array, marray, mvalue)

 其中mvalue原则上为scalar，如果是多维的，则维度需要与array匹配，或者与array的右侧几维匹配。此函数的功能是把array中对应位置的marray不等于mvalue的部分设为缺测值，即mask掉第2和第3个参数不等的部分。这个思路就十分离谱，比如我想要mask掉波作用通量Fx中西风风速uwnd小于5的部分，我需要这么写：

1
Fx = mask(Fx, uwnd.lt.5, False)

 结论就是使用时，第二个参数写上想要mask的条件判断式，第三个参数写False。。。

clmDayTLL/clmDayTLLL, calcDayAnomTLL

 用于计算逐日气候态（年循环）和从原数据中去除年循环。

 前两个函数的主要差别就是一个接收三维数据一个接收思维数据，除原数据之外还需要提供一个yyyyddd数据，对应输入数据的时间维，其中ddd指的是”在一年中是第几天“，取值范围为001-365。计算yyyyddd可以通过cd_calendar(time, -2)接yyyymmdd_to_yyyyddd实现，也可以cd_calendar(time, 0)然后提取结果的(:,0:2)再使用day_of_year(year, month, day)实现。

 calcDayAnomTLL接收一个clmDayTLL函数输出的结果作为参数，操d的是这个函数他…没有四维版本！所以四维数据要转换为三维处理完再reshape回去。多说无益，上一个我实际使用的例子：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
hgt_tmp   = f->hgt;(:,{1000,850,500,200},:,:)
lat       = f->lat
lon       = f->lon
timetmp   = f->time

timehgt     = cd_calendar(timetmp, -2) ;return YYYYMMDD in integer
yyyyddd     = yyyymmdd_to_yyyyddd(timehgt)

dims_hgt = dimsizes(hgt_tmp)
ntime    = dims_hgt(0)
nlev     = dims_hgt(1)
nlat     = dims_hgt(2)
nlon     = dims_hgt(3)

hgt_tmp2 = reshape(hgt_tmp, (/ntime,nlev,nlat*nlon/)) ;TLLL->TLL as required by function
hgt_tmp2!0 = "time"
hgt_tmp2!1 = "level"
hgt_tmp2!2 = "lat_and_lon"
hgt_tmp2&time   = hgt_tmp&time
hgt_tmp2&level  = hgt_tmp&level
printVarSummary(hgt_tmp2)

hgt_clm_366 = clmDayTLL(hgt_tmp2, yyyyddd)
hgt = reshape(calcDayAnomTLL(hgt_tmp2, yyyyddd, hgt_clm_366), dims_hgt)
copy_VarCoords(hgt_tmp, hgt)
printVarSummary(hgt)

此处不得不再次感叹ncl的博大精深！真是太tm有想法了！

填色中央设置为透明

 最后再献丑给大家分享一个自己写的段子，主要功能是把等值线填色图靠近中间值的几个level填色设置成透明，这样就可以省去填色等值分开画还要overlay的麻烦。原理也很简单：把rgba值的最后一个值设为0，因此rgb格式的色表需要手动加上最后一维（值均设为255.0）

 透明半径由set_trans_radius = n指定，如果填色总格数是偶数，则两侧各涂n格；如果为奇数，则中间一格及左右两侧各n-1格设为无色。例如：一共14格填色，n=1，则第7和第8格设为无色；一共15格填色，n=1，则第8格被设为无色。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
res@cnLevelSelectionMode        = "ExplicitLevels"

    lv = fspan(-0.075, 0.075, 11)
    cl = ispan(20, 108, 8)

    colors_tmp = read_colormap_file("MPL_BrBG")
    dim_colors = dimsizes(colors_tmp)

    colors_tmp!0 = "index"
    colors_tmp&index = ispan(0, dim_colors(0)-1, 1)
    colors_tmp!1 = "rgba"
    colors_tmp&rgba = ispan(1, 4, 1)

    colors_work = (/colors_tmp({cl},:)/)

res@cnLevels                    = lv
res@cnFillColors                = colors_work(:,:)

set_trans_radius    = 1         ;设置透明半径

if (dimsizes(cl)%2.eq.0) then
    set_trans = (dimsizes(cl)/2)-1
    set_trans_start = set_trans-set_trans_radius+1
    set_trans_end   = set_trans+set_trans_radius  

    colors_work(set_trans_start:set_trans_end,:) = conform_dims((/set_trans_radius*2,4/), (/255.0,255.0,255.0,0.0/)/255.0, 1)
else
    set_trans = (dimsizes(cl)-1)/2
    set_trans_start = set_trans-set_trans_radius+1
    set_trans_end   = set_trans+set_trans_radius-1  

    colors_work(set_trans_start:set_trans_end,:) = conform_dims((/set_trans_radius*2,4/), (/255.0,255.0,255.0,0.0/)/255.0, 1)
end if
res@cnFillColors                = colors_work(:,:)

 效果图如下：

【物理概念类】

波作用通量（Wave Activity Flux）:

 常用的是Takaya & Nakamura在2001年提出的版本，无需时间平均，公式如下：

 在做合成分析的时候，高度场可以采用合成好的异常场，但是U和V要注意还是使用平均背景场的，是的看走眼用了异常风的就是上周的我！

 顺便推荐一个北海道大学zemi的讲义，总结了各代波作用通量的定义，不过内容是日文，地址点我

垂直水汽积分散度 （VIMFC）

 说到这个又是一个虚惊一场的故事。某日算完该量画完图，突然觉得水汽辐合的位置不对（负的水汽散度代表辐合），检查程序发现计算时竟然乘以了一个-1，一下子头皮发麻…不会是算错了吧！赶紧检查公式：

 正的VIMFC对应着净的降水，即大气失去水分，但是从计算过程看来，这个-1一乘到底还是把辐合辐散给对调了，毕竟垂直积分uq和vq的时候可没有乘以-1.具体情况之后再研究，附上一个示例脚本链接。

不管怎么说，我这一周算是放假啦！