![[PukiWiki]](image/wikitop2.png "[PukiWiki]")
例：流線関数偏差と波活動度フラックス †2021年1月7日～11日の5日平均の250hPa流線関数偏差をTakaya and Nakamura (2001)の波活動度フラックスとともに示す。この期間で平均した東西・南北風速，およびそれらの気候値はすでに計算されているものとする。  スクリプト †begin
;;; データを読む
;;; uv.ncには5日移動平均の250hPa面のu, vおよび対応する気候値uc. vcが入っているとする
in = addfile("uv.nc","r") ; ファイルを開く
lat = in->lat ; 緯度
print(lat(0)) ; 球面調和関数ルーチンを使用するため，
; 緯度座標が南から北であることを確認
u = in->u ; 東西風
uc = in->uc ; 東西風気候値
ua = u-uc ; 東西風偏差
copy_VarCoords(u,ua) ; 東西風偏差に座標情報をコピー
v = in->v ; 南北風
vc = in->vc ; 南北風気候値
va = v-vc ; 南北風偏差
copy_VarCoords(v,va) ; 南北風偏差に座標情報をコピー
dum = uv2sfvpF(ua,va) ; 風偏差から流線関数と速度ポテンシャルを計算
sf = dum(0,:,:) ; 流線関数偏差のみ取り出す
delete(dum) ; 消しておく
;;; Takaya and Nakamura (2001)の波活動度フラックスの定常成分を計算
;;; 微分に使うgradsfが，全球を覆う等間隔格子であること，
;;; 緯度が南から北へ入っていること，
;;; 関数ではなくprocedureであることに注意
dx = sf ; 流線関数のx微分が入る配列を作成
dy = sf ; 流線関数のy微分が入る配列を作成
gradsf(sf,dx,dy) ; 微分（勾配）の計算
dxdx = dx ; dxのx微分が入る配列を作成
dxdy = dx ; dxのy微分が入る配列を作成
gradsf(dx,dxdx,dxdy) ; 微分（勾配）の計算
dydx = dy ; dyのx微分が入る配列を作成
dydy = dy ; dyのy微分が入る配列を作成
gradsf(dy,dydx,dydy) ; 微分（勾配）の計算
d2r = get_d2r("float") ; 度数法から弧度法への変換係数（つまりπ/180）
p = 250./1000. ; 1000hPaで規格化した気圧
; 鉛直座標をlog-p座標で定式化した場合に必要となるスケーリングファクターであって，
; p座標で定式化した場合には表れない
; いまは特定の気圧面でしか描かないので必須ではないと思う
coef = 0.5*p*conform_dims(dimsizes(u),cos(lat*d2r),0)/sqrt(uc*uc+vc*vc)
; 波活動度フラックス全体にかかる係数
wafx = coef*(uc*(dx*dx-sf*dxdx)+vc*(dx*dy-sf*dxdy)) ; 波活動度フラックスのx成分
wafy = coef*(uc*(dx*dy-sf*dxdy)+vc*(dy*dy-sf*dydy)) ; 波活動度フラックスのy成分
sf = sf*1.e-6 ; プロットのために流線関数偏差を10^6でスケーリング
copy_VarCoords(u,wafx) ; 座標情報のコピー
copy_VarCoords(u,wafy)
copy_VarCoords(u,sf)
;;; ここからお絵描きの設定
res1 = True ; 気候値東西風用のres
res1@gsnDraw = False ; gsn_csm_contour実行時にDrawしない
res1@gsnFrame = False ; gsn_csm_contour実行時にFrame更新しない
res1@cnFillOn = True ; シェードを描く
res1@cnLinesOn = False ; 等値線を描かない
res1@cnInfoLabelOn = False ; 等値線の情報を描かない
res1@cnLineLabelsOn = False ; 等値線のラベルを描かない
res1@mpMinLatF = -15 ; 地図の南端の緯度
res1@mpMaxLatF = 75 ; 地図の北端の緯度
res1@mpMinLonF = -90 ; 地図の西端の経度
res1@mpMaxLonF = 270 ; 地図の東端の経度
res1@mpCenterLonF = 90 ; 地図の中央の経度
res1@gsnAddCyclic = True ; 東西方向に周期的
; （これを入れないと東経0度に沿って絵が欠ける）
res1@tmXBLabelFontHeightF = 0.012 ; ヨコ軸ラベルの文字の大きさ
res1@tmYLLabelFontHeightF = 0.012 ; タテ軸ラベルの文字の大きさ
res1@tmBorderThicknessF = 1.25 ; 図の枠線の太さ
res1@tmXBMajorThicknessF = 1.25 ; ヨコ軸の主目盛の太さ
res1@tmYLMajorThicknessF = 1.25 ; タテ軸の主目盛の太さ
res1@pmTickMarkDisplayMode = "Always" ; 軸ラベルをいい感じにするおまじない
res1@mpGeophysicalLineThicknessF = 0.75 ; 地図の線の太さ
res1@mpGeophysicalLineColor = "darkgreen" ; 地図の線の色
res1@pmLabelBarWidthF = 0.62 ; カラーバーの幅
res1@pmLabelBarHeightF = 0.06 ; カラーバーの高さ
res1@pmLabelBarOrthogonalPosF = 0.35 ; カラーバーのタテ方向の位置
res1@pmLabelBarParallelPosF = 0.38 ; カラーバーのヨコ方向の位置
res1@lbLabelFontHeightF = 0.014 ; カラーバーのラベルの文字の大きさ
res1@lbBoxEndCapStyle = "TriangleBothEnds" ; カラーバーの両端を三角形にする
res1@lbTitleOn = True ; カラーバーのタイトルをつける
res1@lbTitleString = "climatological zonal wind [m/s]" ; タイトルの文字列
res1@lbTitleFontHeightF = 0.014 ; カラーバーのタイトルの文字の大きさ
res1@lbBoxLineThicknessF = 1. ; カラーバーの枠線の太さ
res1@gsnLeftString = "2021/01/07-01/11" ; 図の左上の文字
res1@gsnLeftStringFontHeightF = 0.022 ; 図の左上の文字の大きさ
res1@gsnRightString = "250hPa streamfunction anom. [10~S~6 ~N~m~S~2~N~/s]"
; 図の右上の文字
; 流線関数はres2で設定するplot2で描画するが，
; その文字情報はここではres1で設定しておく
res1@gsnRightStringFontHeightF = 0.015 ; 図の右上の文字の大きさ
res1@cnLevelSelectionMode = "ManualLevels" ; 等値線の設定をマニュアルにする
res1@cnLevelSpacingF = 10 ; 等値線（ここではシェード）の間隔
res1@cnMinLevelValF = 0 ; 等値線の最小値
res1@cnMaxLevelValF = 70 ; 等値線の最大値
color = read_colormap_file("sunshine_9lev") ; カラーマップ「sunshine_9lev」をRGBα配列で取得
color(0,:) = (/.600,.600,.600,1.00/) ; 0以下（東風領域）に相当する最も左側の色を灰色（完全不透明）に置き換え
res1@cnFillPalette = color ; カラーマップの設定
res2 = True ; 流線関数偏差用のres
res2@gsnDraw = False ; gsn_csm_contour実行時にDrawしない
res2@gsnFrame = False ; gsn_csm_contour実行時にFrame更新しない
res2@cnFillOn = False ; シェードを描かない
res2@cnLinesOn = True ; 等値線を描く
res2@cnInfoLabelOn = False ; 等値線の情報を描かない
res2@cnLineLabelsOn = True ; 等値線のラベルを描く
res2@cnLineLabelPlacementMode = "constant" ; 等値線ラベルのアルゴリズムを"constant"に設定
res2@cnLineDashSegLenF = 0.10 ; 等値線の破線のセグメントの長さを変える
; 等値線ラベルの間隔を調節するために変える
; cnLineLabelPlacementMode = "constant"のときに有効な手法
res2@cnLineLabelFontHeightF = 0.008 ; 等値線のラベルの文字の大きさ
res2@cnLineThicknessF = 1. ; 等値線の太さ
res2@gsnContourNegLineDashPattern = 11 ; 負の等値線の線種（破線）
; 上で変更したcnLineDashSegLenFに応じて破線の間隔が変わってしまうことに注意
res2@gsnContourZeroLineThicknessF = 0 ; 0の等値線の太さ（0なので描かない）
res2@gsnContourPosLineDashPattern = 0 ; 正の等値線の線種（実線）
res2@gsnAddCyclic = True ; 東西方向に周期的
; （これを入れないと東経0度に沿って絵が欠ける）
res2@cnLevelSelectionMode = "ManualLevels" ; 等値線の設定をマニュアルにする
res2@cnLevelSpacingF = 6 ; 等値線の間隔
res2@cnMinLevelValF = -60 ; 等値線の最小値
res2@cnMaxLevelValF = 60 ; 等値線の最大値
res2@cnLineDrawOrder = "PostDraw" ; 通常のdrawの後から描く（地図の線よりも手前に出すため）
res3 = True ; 波活動度フラックス用のres
res3@gsnDraw = False ; gsn_csm_vector実行時にDrawしない
res3@gsnFrame = False ; gsn_csm_vector実行時にFrame更新しない
res3@vcGlyphStyle = "FillArrow" ; 矢印の描画スタイルとしてFillArrowを用いる
res3@vcFillArrowFillColor = "blue" ; 矢印の塗りつぶしの色
res3@vcFillArrowWidthF = 0.04 ; 矢印の太さ
res3@vcFillArrowHeadXF = 0.4 ; 矢印の頭の長さ（外側）
res3@vcFillArrowHeadYF = 0.12 ; 矢印の頭の幅
res3@vcFillArrowHeadInteriorXF = 0.2 ; 矢印の頭の長さ（内側）
res3@vcRefLengthF = 0.04 ; 基準となる矢印の長さ
res3@vcRefMagnitudeF = 75 ; 基準となる矢印に対応する値
wafx = where(sqrt(wafx*wafx+wafy*wafy).lt.res3@vcRefMagnitudeF/5, \
wafx@_FillValue,wafx)
; 基準となる矢印の1/5に満たない長さのベクトルは未定義値にしてしまい描かない
; 未定義値を入れるのはwafxかwafyのどちらでもよい
wafy = where(uc.lt.0,wafy@_FillValue,wafy)
; 気候値東風領域には描かない
; 未定義値を入れるのはwafxかwafyのどちらでもよい
res3@vcRefAnnoOn = True ; 矢印の凡例をつける
res3@vcRefAnnoFontHeightF = 0.010 ; 矢印の凡例の文字の大きさ
res3@vcRefAnnoString1 = "wave activity flux" ; 凡例上側の文字列
res3@vcRefAnnoString2On = True ; 矢印の凡例下側に文字を書く
res3@vcRefAnnoString2 = res3@vcRefMagnitudeF+" [m~S~2~N~/s~S~2~N~]"
; 凡例下側の文字列
res3@vcRefAnnoPerimOn = False ; 凡例に枠線をつけない
res3@vcRefAnnoOrthogonalPosF = -0.38 ; 矢印の凡例のタテ方向の位置
res3@vcMinDistanceF = 0.024 ; 矢印の間隔の最小値
res3@vcVectorDrawOrder = "PostDraw" ; 通常のdrawの後から描く（地図の線よりも手前に出すため）
wks = gsn_open_wks("eps","sfwaf") ; workstation の定義
plot1 = gsn_csm_contour_map_ce(wks,uc,res1) ; 気候値東西風の絵を作成
plot2 = gsn_csm_contour(wks,sf,res2) ; 流線関数偏差の絵を作成
plot3 = gsn_csm_vector(wks,wafx,wafy,res3) ; 波活動度フラックスの絵を作成
overlay(plot1,plot2) ; 重ねる
overlay(plot1,plot3) ; 重ねる
draw(plot1) ; 描画を実行
frame(wks) ; wksの更新
end
|
