NCLでキレイな絵が描けるようになると,アニメーションを作りたくなってくる。しかし,NCLには直接的にアニメーションを作成する機能はない。したがって,まずNCLで絵をたくさん描いてから,それらを別のツールで一つのアニメーションに変換することが必要となる。
それなりのアニメーションを作るには,数十~数百枚の絵を用意しなければならないと思うが,この量の絵を普通に生成するのにはそれなりの時間がかかることもあるだろう。そこで,ここでは,NCLスクリプト内で別の複数のサブプロセスを並列実行できる関数であるsubprocessを活用して,効率的に絵を生成し,最後にimagemagickを利用してgifアニメーションを作成する例を紹介する。
※クリックでgifアニメーションのページに移ります
JRA-55の6時間値(GRIB形式)から2019年1月の850hPa気温のアニメーションを作成する。NCLスクリプトは,データを読んで絵を描くdraw.nclと,draw.nclを並列実行するdriver.nclに分かれる。
変数stepおよびddhhは,実行時にdriver.nclから受け取る。
begin dir = "./work/" ; 絵を吐き出すディレクトリ infile = "/data/JRA-55/Hist/Daily/anl_p125/201901/anl_p125_tmp.201901"+ddhh ; 読むデータ,ddhhは日時が入った文字列。 in = addfile(infile+".grb","r") ; ファイルを開く。拡張子を付与。 tmp = in->TMP_GDS0_ISBL({850},:,:) ; 850hPa面の気温を読む ;; pngで図を作成,ファイル名はfig_???.pngのように時間ステップ順に連番とする wks = gsn_open_wks("png",dir+"fig_"+sprinti("%0.3i",step)) res = True res@cnFillOn = True ; シェードを描く res@cnLinesOn = True ; 等値線を描く res@cnInfoLabelOn = False ; 等値線の情報なし res@cnLineLabelsOn = False ; 等値線のラベルなし res@cnLineThicknessF = 0.8 ; 等値線の太さ res@cnLineColor = "white" ; 等値線の色 res@cnFillPalette = "NCV_bright" ; シェードに使うカラーパレット res@pmLabelBarWidthF = 0.85 ; カラーバーの幅 res@lbLabelFontHeightF = 0.016 ; カラーバーのラベルの文字の大きさ res@pmLabelBarOrthogonalPosF = 0.2 ; カラーバーのY方向の位置 res@lbBoxEndCapStyle = "TriangleBothEnds" ; カラーバーの端を両側三角にする res@lbTitleOn = True ; カラーバーのタイトルを書く res@lbTitleString = tmp@units ; カラーバーのタイトル(ここでは気温の単位) res@lbTitleFontHeightF = 0.018 ; カラーバーのタイトルの文字の大きさ res@lbTitlePosition = "Right" ; カラーバーのタイトルの位置 res@lbTitleDirection = "Across" ; カラーバーのタイトルの文字の向き res@tmBorderThicknessF = 1.5 ; 図の枠の線の太さ res@tmXBMajorThicknessF = 1.5 ; X軸の主目盛の線の太さ res@tmYLMajorThicknessF = 1.5 ; Y軸の主目盛の線の太さ res@mpGeophysicalLineColor = "gray20" ; 地図の線の色 res@mpGeophysicalLineThicknessF = 1.5 ; 地図の線の太さ res@pmTickMarkDisplayMode = "Always" ; 軸ラベルをいい感じにする res@tmXBLabelFontHeightF = 0.018 ; X軸のラベルの文字の大きさ res@tmYLLabelFontHeightF = 0.018 ; Y軸のラベルの文字の大きさ res@gsnLeftStringFontHeightF = 0.03 ; 図の左上の文字の大きさ res@gsnLeftString = "850hPa Temp. (201901"+ddhh+")" ; 図の左上の文字列 res@gsnRightString = "" ; 図の右上の文字列(なし) res@cnLevelSelectionMode = "ManualLevels" ; 等値線間隔をマニュアルモードに res@cnLevelSpacingF = 3 ; 等値線間隔 res@cnMaxLevelValF = 303 ; 等値線最大値 res@cnMinLevelValF = 243 ; 等値線最小値 plot = gsn_csm_contour_map_ce(wks,tmp,res) ; 正距円筒図法で描画 end
draw.nclを適切に並列実行させて効率的に絵を生成し,最後にトリミングとアニメーションの作成を行う。
;;; ;;; サブプロセスを並列実行するprocedureであるpara_processをはじめに定義する。 ;;; もちろん別スクリプトに記述して,はじめにloadしてもよい。 ;;; 実行するコマンド群を文字列配列commandとして,最大並列実行数を整数MAX_CONCURRENTとして受け取る。 ;;; MAX_CONCURRENTは計算機のコア数や(共用計算機なら)他人への迷惑を考慮しつつ決めましょう。 ;;; undef("para_process") procedure para_process(command[*]:string,MAX_CONCURRENT[1]:integer) local numCompleted, numActive, numSteps, pid, i ;; 内部変数 begin numSteps = dimsizes(command) ;; 実行するプロセスの総数 numCompleted = 0 ;; 終了したプロセス数 numActive = 0 ;; 実行中のプロセス数 do i = 0, numSteps-1 do while(numActive.ge.MAX_CONCURRENT) ;; 実行中のプロセス数が与えた最大実行数に達している間は以下の処理を行う pid = subprocess_wait(0,True) ;; いずれかのプロセスが終了するまで待つ if (pid.gt.0) then ;; もしあるプロセスが終了したなら... numCompleted = numCompleted + 1 ;; 完了したプロセス数を1増やし numActive = numActive - 1 ;; 実行中のプロセス数を1減らし break ;; do whileを抜ける end if end do pid = subprocess(command(i)) ;; i番目のコマンドを実行 numActive = numActive + 1 ;; 実行中のプロセス数を1増やす end do ;; ここに来た頃にはすべてのプロセスを実行済(実行中か完了後)である do while (numCompleted.lt.numSteps-1) ;; まだすべてのプロセスが完了していない間は以下の処理を行う pid = subprocess_wait(0,True) ;; いずれかのプロセスが終了するまで待つ if (pid.gt.0) then ;; もしあるプロセスが終了したなら完了したプロセス数を1増やす numCompleted = numCompleted + 1 end if end do ;; すべてのプロセスが完了 end ;;; ここからがメイン begin MAX_CONCURRENT = 10 ;; サブプロセスの最大並列実行数 yyyymmddhh = yyyymmddhh_time(2019,2019,6,"integer") ;; 2019年のYYYYMMDDHHを整数で取得 JAN = ind(yyyymmddhh.lt.2019020000) ;; YYYYMMDDHHが2019020000より小さい,つまり1月の部分のみの部分に対応する配列 ddhh = sprinti("%0.4i",yyyymmddhh(JAN)%10000) ;; yyyymmddhhから1月の部分のみ切り出し,10000で割った余りを計算することでDDHHを得る ;; ついでに文字列に変換(3桁の場合は0で埋める) ddhh = str_get_dq()+ddhh+str_get_dq() ;; ダブルクォートで括る。ただし,スクリプト上で " を文字列として ;; 用いることができないので関数str_get_dqを用いている。 ;; つまりddhhは「"1206"」のような文字列 numSteps = dimsizes(ddhh) ;; ddhhの長さ,サブプロセスの総数となる nclcomm = "ncl -nQ draw.ncl" ;; draw.nclを実行するコマンド(の一部) trmcomm = "convert -trim +repage" ;; 画像をトリミングするコマンド(の一部) last = sprinti("%0.3i",numSteps-1) gifcomm = "convert -layers optimize -loop 0 -delay 25 ./work/fig_*.png -delay 75 ./work/fig_"+last+".png anim.gif" ;; 連番pngファイルからgifアニメーションを作成するコマンド ;; -loop 0 はループ再生を,-delay 25 は0.25秒間隔を意味する ;; 最後の絵だけ,-delay 75 で0.75秒追加している I = ispan(0,numSteps-1,1) ;; 0からnumSteps-1までの整数の配列を作成 ;;; 上で定義したpara_processを通して,numSteps個のdraw.nclを実行 ;;; これは ncl -nQ draw.ncl 'step=13' 'ddhh="0406"' のようなコマンドである para_process(nclcomm+" 'step="+I+"' 'ddhh="+ddhh+"'",MAX_CONCURRENT) ;;; 絵の生成がすべて完了したら, ;;; 上で定義したpara_processを通して,numSteps個のトリミングコマンドを実行 png = "./work/fig_"+sprinti("%0.3i",I)+".png" ;; pngファイルの配列 para_process(trmcomm+" "+png+" "+png,MAX_CONCURRENT) ;;; 絵のトリミングがすべて完了したら,gifアニメーションを作成するコマンドを実行 system(gifcomm) end