LoginSignup
4
1
@e_toyoda(英司 豊田)

TFP (thermal front parameter) で前線ぽいものを描画する

Last updated at Posted at 2023-11-23

数値予報GPVを描画して遊んでいるわけですが、前線などのオブジェクトも自動判定して描いてみたいと思いました。

定義1

実際に業務利用されて明瞭にアルゴリズムが書かれている例として、[1] では次の三条件を満たすものをNP前線と呼んでいます。

\nabla_n {\rm DDT} = 0 \quad\cdots (1)

{\rm DDT} \ge 0.5 {\rm K} / \left(381 {\rm km}\right)^2 \quad\cdots (2)

\nabla_n\cdot\nabla_n {\rm DDT} \le 0 \quad\cdots (3)

ここで $T$ は850 hPa面気温、$\nabla_n$ は $T$ の増大する方向の微分、

{\rm DDT} = -\nabla_n\cdot \left| \nabla_n\cdot T \right|
\quad\cdots(4)

であって、温暖前線と寒冷前線の区別は1000 hPa面地衡流の向きや大きさによるのだといいます。

[1] 予報部数値予報課(1988): 会話型気象図作成支援処理システムについて. 測候時報55.4 pp.181-191.

定義2

もうちょっと歴史を通覧したレビューが [2] にあって、Renard and Clarke (1965) [3] の thermal front parameter (TFP) というのが始祖になっていて、そこから色々チューニングしているらしいということがわかります。

[2] 北畠尚子(2023): 総観気象学 応用編. 気象庁. https://www.jma.go.jp/jma/kishou/know/expert/
[3] Renard, R. and L. Clarke (1965): Experiments in numerical objective frontal analysis. Mon Wea. Rev., 93, 547-556. https://doi.org/10.1175/1520-0493(1965)093%3C0547:EINOFA%3E2.3.CO;2

今日TFPと呼ばれるものは、なんらかの熱的変数（気温、相当温位など）を $\tau$ として

{\rm TFP} = -\nabla \left|\nabla\tau\right|\cdot \frac{\nabla\tau}{\left|\nabla\tau\right|}
\quad\cdots(5)

として算出したもので、ざっくりこれを極大化するところを前線とするのです。ただ極大だけだと要らんものがたくさんひっかかるので、定義1の三条件 (1)...(3) のような選別をします。アイデアは、温度傾度集中帯の南端を前線とする伝統的なスタイルを再現するため、TFP をさらに微分してゼロ点を探し（典型的な前線構造では前線帯の前後に2つ現れる）、さらに微分して負値となる条件を課して暖気側を選ぶということになります。

いっぽうで、近年は高解像度の分析では高次の微分によるノイズを防ぐため、 ${\rm TFP} = 0$ を用いることもあるそうです [2]。その場合は (2) 式に代えて、一次微分の絶対値 $\left|\nabla\tau\right|$ による足切りなどが必要になります。

もっと簡単化できんのか

ところでヘンテコな微分がたくさん出てきましたが、要は前線断面に沿って2回微分しているだけです。それならいっそのこと単なるラプラシアン

{\rm div}\ {\rm grad}\ \tau = \nabla^2\tau = 0
\quad\cdots (6)

なる線を描いてみたらどうなんだろうという気がしてきますね。

実際にやってみる

コード置場など

で、実際にやってみたので、quick hack code の説明を置いておきます。コードはこのへんです
https://github.com/etoyoda/grib2png/issues/36

このプログラムは普段、GISC Tokyo ウェブサイトからダウンロードした 1.25°解像度 GSM 予報値を使ってウェブメルカトル図法の全球 1024x1024 ピクセル画像を作っています。
微分演算はウェブメルカトル格子で行うことにしました。そうしたかった直接の動機は、プログラムの作りから、TFP 計算サブルーチンを生やす場所がそこがやりやすかった、という偶有的な理由なのですが、 メルカトル図法は等角図法なのでベクトル解析が差分で綺麗に書ける という大切な意味もあります。MSM などの領域モデルがランベルト正角図法を使うのと同じ理由ですね。

問題のTFP算出はこのへんです。
https://github.com/etoyoda/grib2png/blame/b3a8d3cbcfb629fc5e3f51c84a896f51a78bc908/visual.c#L546

visual.c:546
  int
drawfront(png_bytep *ovector, const double *gbuf,
  size_t owidth, size_t oheight, palette_t pal)
{

ウェブメルカトル座標の格子に入った浮動小数点値 gbuf (owidth×oheight 個要素の配列) から前線ぽいものの位置を判定して、libpngが出力するラスタデータovector（こちらは4バイトがowidth個並んだ配列へのポインタがoheight個ならんだベクタ）の該当ピクセルに書き込む、ということをします。

先行スムージング

変数宣言とかメモリ確保とかは略。まずはgbufに7x7と5x5の二回移動平均をかけます。入力gbufはより粗い格子系（最も密な熱帯でgbufの1/3程度の格子間隔）から線形補間で作っているので、もともと空間方向の傾度 ${\rm grad}\ \theta_{\rm e}$ が階段関数になっているのですね。これで安直に二回微分をかけるとガタガタになってしまうのでしかたがないのです。

visual.c:562
  for (size_t j = 3; j < oheight - 3; j++) {
    for (size_t i = 0; i < owidth; i++) {
      double sum;
      sum = 0.0;
      for (size_t jc = j - 3; jc <= j + 3; jc++) {
        for (int ics = i - 3; ics <= i + 3; ics++) {
          sum += gbuf[(ics % owidth)+jc*owidth];
        }
      }
      xbuf[i+j*owidth] = sum / 49.0;
    }
  }

これは7x7格子の移動平均です。見ての通り東西はサイクリックに参照できますが南北は幅3格子ぶんだけ算出をあきらめています。まあ所詮移動平均なので元データをコピーしてもだいたい同じです（こういうところがクイックハック）。

傾度絶対値の算出（準備）

visual.c:586
  // dgbuf に傾度ベクトル長を設定
  for (size_t j = 1; j < oheight - 1; j++) {
    for (size_t i = 0; i < owidth; i++) {
      size_t ip1 = (i + 1) % owidth;
      size_t im1 = (i - 1) % owidth;
      dgbuf[i + j*owidth] = 
        hypot(sbuf[ip1+j*owidth] - sbuf[im1+j*owidth],
        sbuf[i+(j+1)*owidth] - sbuf[i+(j-1)*owidth]);
    }
  }

上は傾度ベクトルの絶対値を算出しています。

{\tt dgbuf} =
\left|\nabla\tau\right|
= \left|\left(
\begin{array}
\ \partial_x \\
\partial_y \\
\end{array}
\right)\tau\right|  
= \sqrt{
\left(\partial_x\tau\right)^2
+ \left(\partial_y\tau\right)^2
}

ということですね

TFPの算出

visual.c:596
  // TFP (傾度ベクトルの方向での傾度の微分) を計算
  for (size_t j = 1; j < oheight - 1; j++) {
    for (size_t i = 0; i < owidth; i++) {
      size_t ip1 = (i + 1) % owidth;
      size_t im1 = (i - 1) % owidth;
      double nx, ny;
      nx = (sbuf[ip1+j*owidth]-sbuf[im1+j*owidth])/dgbuf[i+j*owidth];
      ny = (sbuf[i+(j+1)*owidth]-sbuf[i+(j-1)*owidth])/dgbuf[i+j*owidth];
      tfpbuf[i + j*owidth] = 
      - nx * (dgbuf[ip1+j*owidth] - dgbuf[im1+j*owidth])
      - ny * (dgbuf[i+(j+1)*owidth] - dgbuf[i+(j-1)*owidth]);
      // ついでにラプラシアン
      xbuf[i+j*owidth] =
        sbuf[ip1+j*owidth] + sbuf[im1+j*owidth]
      + sbuf[i+(j+1)*owidth] + sbuf[i+(j-1)*owidth] - 4*sbuf[i+j*owidth];
    }
  }

ここがTFPですが、(5) を

{\rm TFP} = - \left(\nabla \left| \nabla\tau \right|\right)\cdot \frac{\nabla\tau}{\left| \nabla\tau \right|}

と解釈して

{\rm TFP} = - \frac{\partial_x\tau}{{\tt dgbuf}}\partial_x{\tt dgbuf}
- \frac{\partial_y\tau}{{\tt dgbuf}}\partial_y{\tt dgbuf}

としているわけですね。ラプラシアン (6) もついでに計算していますが配列をレジスタみたいに使いまわしています。これを動かしているサーバがメモリ足りないので。

さて、格子の差分を格子間隔で割るのが本当で、メルカトル図法なので緯度によってちょっと（大きく）異なるのですが、そこは端折っています。だから quick hack ですね。

ゼロ線探索と描画

さてゼロ線を探索するんですが、真面目に線を探索しないで、結果の画像に線が描かれさえすればよいというスタンスで、隣接ピクセルと符号が違うピクセルを選んで色を塗っています。
どんだけの太さの線をひきたいかと、千切れを防ぎたいかのかねあいで、8隣接のかわりに4隣接にしたり、ここでは色を塗るセルの TFP を正の値に限定していますが、正負を許して太線にするなどのこともできます。

visual.c:621-658
  for (size_t j = 1; j < oheight - 1; j++) {
    for (size_t i = 0; i < owidth; i++) {
      size_t ip1 = (i + 1) % owidth;
      size_t im1 = (i - 1) % owidth;
      png_bytep pixel = ovector[j] + i * 4;
      if (
        (dgbuf[i+j*owidth] > mingrad) &&
        (tfpbuf[i+j*owidth] > 0.0) &&
        (
          (tfpbuf[im1+(j-1)*owidth] * tfpbuf[i+j*owidth] < 0.0) ||
          (tfpbuf[i  +(j-1)*owidth] * tfpbuf[i+j*owidth] < 0.0) ||
          (tfpbuf[ip1+(j-1)*owidth] * tfpbuf[i+j*owidth] < 0.0) ||
          (tfpbuf[im1+j    *owidth] * tfpbuf[i+j*owidth] < 0.0) ||
          (tfpbuf[ip1+j    *owidth] * tfpbuf[i+j*owidth] < 0.0) ||
          (tfpbuf[im1+(j+1)*owidth] * tfpbuf[i+j*owidth] < 0.0) ||
          (tfpbuf[i  +(j+1)*owidth] * tfpbuf[i+j*owidth] < 0.0) ||
          (tfpbuf[ip1+(j+1)*owidth] * tfpbuf[i+j*owidth] < 0.0)
        )
      ){
        pixel[0] = 128; pixel[1] = pixel[2] = 12; pixel[3] = 255;
      }
    }
  }

足切りの mingrad ですが、見た感じそれっぽくなるように適当にチューニングして、さしあたり次のようにしています。

  • 気温の場合 12.0
  • 相当温位の場合 32.0

これは、2ピクセルあたりの温度差（処理の都合上0.1K単位）ですから、たとえば北緯30度の場合は $2\cdot 40,000{\rm km}\cdot{\rm cos}30^\circ/1024 = 67.7{\rm km}$ で割って

  • 気温の場合 1.77 K/100km
  • 相当温位の場合 4.73 K/100km

ということになりますね。Jonkner et al [4] は相当温位で 4.5 K/100km を使ったので、まあそんなもんなのでしょう。

[4] Jenkener et al. 2009: Detection and climatology of fronts in a high-resolution model reanalysis over the Alps https://doi.org/10.1002/met.142

ラプラシアンで代用するのはやっぱり駄目

TFPの微分演算をラプラシアンで代用したらどうなるかというのは、試してみて、やっぱりこりゃ駄目だと思いました。前線然とした方向以外の微分を考慮するわけですから、そりゃ蛇行します。
まあいうてそれほどダメでなければ div grad という物理的意味も明確なのでいいかなーと思ったのですが、実際試してみてダメだった、というのが [3] の式を鵜呑みにするまえにやりたかったことです。
https://twitter.com/e_toyoda/status/1728220396468736104

4
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
1