![[PukiWiki]](image/wikitop2.png "[PukiWiki]")
例：ある観測点の気温と降水量の季節変化 †気象庁のページからダウンロードした神奈川県三浦の気温と降水量について，1981年～2010年の気候値を計算し，XYプロットで描画。 
ダウンロードした時刻：2017/12/06 16:18:35
, 三浦, 三浦, 三浦
年月日 , 平均気温(℃), 平均気温(℃), 平均気温(℃)
, , 品質情報, 均質番号
1981/1/1, 5.9, 8, 1
1981/1/2, 6.7, 8, 1
1981/1/3, 5.5, 8, 1
1981/1/4, 6.4, 8, 1
1981/1/5, 4.8, 8, 1
…
スクリプト †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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
out = "eps"
y1 = 1981 ; 最初の年
y2 = 2010 ; 終わりの年
NY = y2-y1+1
YYYYMMDD = yyyymmdd_time(y1,y2,"integer") ; y1からy2までの日付をYYYYMMDDの形で取得
YYYYMM = YYYYMMDD/100 ; YYYYMMの配列を取得
YYYY = YYYYMM/100 ; 年(YYYY)の配列を取得
MM = YYYYMM - YYYY*100 ; 月(MM)の配列を取得
DD = YYYYMMDD - YYYYMM*100 ; 日(DD)の配列を取得
;;; データを読む。ここでは気象庁からダウンロードしたcsvファイルを読むことにする
;;; 日平均気温のデータを読む
data = asciiread("data_temp.csv",-1,"string") ; csvファイルの中身を文字列として取得。行数の大きさを持つ1次元配列
data := data(5:) ; 最初の5行を捨てる
data := str_split_csv(data,",",0) ; ,で区切る
temp0 = tofloat(data(:,1)) ; 各行の2列目を取り出してfloat型に変換
delete([/data/])
;;; 日最高気温のデータを読む
data = asciiread("data_temp_max.csv",-1,"string")
data := data(5:)
data := str_split_csv(data,",",0)
tmax0 = tofloat(data(:,1))
delete([/data/])
;;; 日最低気温のデータを読む
data = asciiread("data_temp_min.csv",-1,"string")
data := data(5:)
data := str_split_csv(data,",",0)
tmin0 = tofloat(data(:,1))
delete([/data/])
;;; 日降水量のデータを読む
data = asciiread("data_prec.csv",-1,"string")
data := data(5:)
data := str_split_csv(data,",",0)
prec0 = tofloat(data(:,1))
delete([/data/])
;;; 時間の次元を年と日付の2次元に変形
DDD = day_of_year(YYYY,MM,DD) ; 年初からの日数をDDDの形で取得
temp = new((/NY,365/),float)
tmax = new((/NY,365/),float)
tmin = new((/NY,365/),float)
prec = new((/NY,365/),float)
do i = 0, NY-1
I = ind(YYYY.eq.i+y1.and.DDD.eq.1) ; 各年について1月1日の場所を取得
temp(i,:) = temp0(I:I+364) ; 各年の1日目から365日目までのデータを切り出し，新しい配列に代入
tmax(i,:) = tmax0(I:I+364) ; なお，閏年の場合は12月31日の値が捨てられていることになる
; (ここでは移動平均するので良しとする)
tmin(i,:) = tmin0(I:I+364)
prec(i,:) = prec0(I:I+364)
end do
delete([/temp0,tmax0,tmin0,prec0/])
temp_clm = dim_avg_n(temp,0) ; 日別気候値を計算
tmax_clm = dim_avg_n(tmax,0)
tmin_clm = dim_avg_n(tmin,0)
prec_clm = dim_avg_n(prec,0)
temp_clm_ann = avg(temp_clm) ; 年平均気候値を計算
tmax_clm_ann = avg(tmax_clm)
tmin_clm_ann = avg(tmin_clm)
prec_clm_ann = avg(prec_clm)
temp_clm = runave(temp_clm,31,-1) ; 31日移動平均を計算する。周期境界条件を用いることができる。
tmax_clm = runave(tmax_clm,31,-1)
tmin_clm = runave(tmin_clm,31,-1)
prec_clm = runave(prec_clm,31,-1)
prec_clm = runave(prec_clm,11,-1) ; 降水量についてはもう一度11日移動平均をかけてなめらかに
;;; ここからは描画設定です
wks = gsn_open_wks(out,"miura")
res1 = True
res1@gsnDraw = False
res1@gsnFrame = False
res1@vpHeightF = 0.33 ; 図の高さを設定
res1@tmXBMode = "Explicit" ; 下部のx軸を"Explicit"モードに変更
;;; "Explicit"モードでは次の二つのresourceでラベルの文字とそれを置く場所を配列として指定する
res1@tmXBLabels = sprinti("%1.0i",ispan(1,12,1)) ; ラベルの文字。1から12までの整数
res1@tmXBValues = DDD(ind(YYYY.eq.2001.and.DD.eq.1))-1
; 毎月1日の年初からの日数を取得し，1を引く
; (1月1日はDDD=1であるが,変数の0番目の値と対応するから)
; ここでは2001年のDDDの値を取っているが，閏年でなければなんでもいい
res1@trXMinF = 0 ; x軸は0～364
res1@trXMaxF = 364
res1@tmXTOn = False ; 上部のx軸は使わない
res1@tmXBLabelFontHeightF = 0.022 ; x軸のラベルの文字の大きさ
res1@tmXBMajorThicknessF = 4 ; x軸の目盛りの太さを太くする
;;; res2を作成。これ以降，res1に降水量のプロット(左側のy軸を用いる)，
;;; res2に気温のプロット(右側のy軸を用いる)の設定を入れていく
res2 = res1
res1@trYMinF = 0 ; 左側のy軸の範囲は0～8(mm/day)
res1@trYMaxF = 8
res1@tmYLMode = "Manual" ; 左側のy軸を"Manual"モードに変更
res1@tmYLMinorPerMajor = 10 ; 大目盛の間に10本の小目盛を入れる
res1@tmYLMinorLineColor = "blue" ; 小目盛の色を青にする
res1@tiYAxisString = "precipitation [ mm/day ]" ; 軸のタイトル
res1@tiYAxisFontHeightF = 0.022 ; 軸のタイトルの文字の大きさ
res1@tiYAxisFontColor = "green4" ; 軸のタイトルの文字の色
res1@tmYLLabelFontHeightF = 0.022 ; ラベルの文字の大きさ
res1@xyLineThicknessF = 3. ; 折れ線の太さ
res1@xyLineColor = "blue" ; 折れ線の色
res1@gsnYRefLine = 0 ; y=0に線を引く(x軸と被っているので見えない)
res1@gsnAboveYRefLineColor = "aquamarine" ; 上で設定した線の上部を青緑で塗りつぶす
res2@trYMinF = 0 ; 右側のy軸の範囲は0～8(mm/day)
res2@trYMaxF = 32
res2@tmYRMode = "Manual" ; 右側のy軸を"Manual"モードに変更
res2@tmYRMinorOn = False ; 小目盛を用いない
res2@tmBorderLineColor = "red" ; 図の枠の色を赤に変更
res2@tmYLBorderOn = False ; 左側の枠線を消す
res2@tmXTBorderOn = False ; 上側の枠線を消す
res2@tmXBBorderOn = False ; 下側の枠線を消す
; これによって右側のy軸部分の枠線のみが赤で描かれ，
; あとでres1で描かれたデフォルトの枠線の上に重ねられることになる
res2@tmYRLabelFontColor = "purple" ; ラベルの文字の色
res2@tmYRLabelFontHeightF = 0.022 ; ラベルの文字の大きさ
res2@tmYRMajorLineColor = "orange" ; 大目盛の色をオレンジに
res2@tiYAxisString = "temperature [ ~S~o~N~C ]" ; 軸のタイトル
res2@tiYAxisFontHeightF = 0.022 ; 軸のタイトルの文字の大きさ
res2@tiYAxisFontColor = "magenta" ; 軸のタイトルの文字の色
res2@xyLineThicknessF = 3 ; 折れ線の太さ
res2@xyLineColor = "red" ; 折れ線の色
;;; 気温と降水量の年平均値を水平な線で示す
;;; しかし，気温のプロットに用いるres2で両方の値を描画することから，
;;; 降水量については気温の軸での値に変換する必要があるので変換する。
prec_clm_ann = res2@trYMinF + (res2@trYMaxF-res2@trYMinF) * \
(prec_clm_ann-res1@trYMinF)/(res1@trYMaxF-res1@trYMinF)
res2@gsnYRefLine = (/temp_clm_ann,prec_clm_ann/) ; それぞれの年平均値のところに水平な線を引く
res2@gsnYRefLineDashPatterns = (/3,5/) ; 年平均値の線種
res2@gsnYRefLineThicknesses = (/2,2/) ; 年平均値の太さ
res2@gsnYRefLineColors = (/"red","blue"/) ; 年平均値の色，気温が赤で降水量が青
res1@gsnLeftString = "Miura, Kanagawa (1981-2010)" ; 図のタイトル(ここではgsnLeftStringで設定)
res1@gsnLeftStringFontHeightF = 0.03 ; タイトルの文字の大きさ
;;; ここでやっとプロット。2つのy軸を用いるgsn_csm_xy2を使用
plot1 = gsn_csm_xy2(wks,ispan(0,364,1),prec_clm,temp_clm,res1,res2)
;;; 最高気温と最低気温を平均気温のまわりに帯を描くことで表現したい
;;; 最高・最低気温の折れ線を左右端でつなげた多角形の内部を塗りつぶした図形を描く
gsres = True
gsres@gsFillColor = "pink" ; シェードの色
gsres@gsFillOpacityF = 0.6 ; 不透明度(0が完全に透明)
gsres@tfPolyDrawOrder = "PreDraw" ; 標準的にDrawで描かれるものより前に描く(平均気温の赤線の背後に位置させたい)
XP = new(365*2,integer) ; x座標が入る配列
YP = new(365*2,float) ; y座標が入る配列
do k = 0, 364
XP(k) = k ; 前半は1/1～12/31に対応する座標
XP(365*2-1-k) = k ; 後半は上の逆順で値を代入
YP(k) = tmax_clm(k) ; 前半は最高気温の値(1/1～12/31の順)
YP(365*2-1-k) = tmin_clm(k) ; 後半は最低気温の値(12/31～1/1の順)
end do
;;; 図形の座標は降水量のy軸の値が基準となる(res1に対応する)ため，気温の値を降水量の軸に対応する値に変換する。
YP = res1@trYMinF + (res1@trYMaxF-res1@trYMinF) * \
(YP-res2@trYMinF)/(res2@trYMaxF-res2@trYMinF)
;;; 図形をplot1上に描画
dum0 = gsn_add_polygon(wks,plot1,XP,YP,gsres)
draw(plot1) ; plot1を描画
frame(wks) ; wksを更新
end
|
