[2023-06-07] 格子点データをめぐって

空間分布する量の表現

• 観測点の位置での数量
• 規則的に分布する格子点の位置での数量
• 関数の重みつき重ね合わせ (「直交関数展開」「スペクトル法」)
  • 平面ならば、直交する2軸それぞれ、三角関数 (sin, cos) によるフーリエ変換
  • 球面では、球面調和関数 (「球関数」)、東西方向は三角関数、南北方向はルジャンドル陪関数

格子点データのはじまり

• 気象では1950年代からシミュレーションのためのモデルに必要とされる
  • 格子点差分モデル
  • スペクトルモデル ... 放射や雲の計算のために格子点も併用
• 等値線表現の中間段階としての格子点データ
• 行政的あるいは地理学的な統計解析のための (日本語でいう) メッシュデータ (cf. ポリゴンデータ)

格子点 (grid point) の値か、ます目 (grid box) の値か? (grid registration vs. pixel registration)

地球上の位置の表現

(球面の表示装置もできてはきたが) 通常は、平面の紙面や画面に表示するので、地図投影が必要。

格子データをどのような座標軸のもとでつくるかと、表示の座標軸とは、かならずしも同じではない。

わたしの教材はいまのところ、格子データも表示も、経度を横軸、緯度を縦軸にとったもの。 地図投影としては「正距円筒図法」、俗に「方眼図法」。英語では Plate Carree (フランス語起源)。

その他の地図投影

• 北極・南極をふくむ地域をあつかうには、極を中心とする方位図法。
• 天気図でよくつかわれるのは、等角投影 ... ローカルな角度が正しい。 (グローバルな方位角が正しいわけではない。)
  • 極中心ならばステレオ図法。
  • 緯線経線を直線とする円筒図法ならばメルカトル図法。(縮尺が緯度によってちがう。熱帯だけあるいは小地域用ならばよい。)
  • 中緯度ではランベルト等角円錐図法がつかわれることがある。
    • ふつうは地軸を軸とするが、試作段階の解析雨量では斜軸がつかわれていた (日本列島にそった領域を分割したかったので)。

[地図投影 (map projection) (地理情報学教材 2000年]

鉛直は、気象では、p座標がつかわれることが多く、等間隔でないことが多い。

3次元の表現は2次元の表現ほど発達していない。 また、大気は連続体で、雲なども明確な 3次元立体の表面をもたないので、 いわゆる3次元グラフィックスよりも、平面をふくむ 2次元断面 で見たほうがよいことが多い。

気象モデルの格子

差分モデルでは、その差分計算方式の必要性できめられる。たとえば風の東西成分と南北成分の位置をたがいちがいにとる。

スペクトルモデルでは、ルジャンドル陪関数の値がゼロになる緯度 (ガウス緯度) にとる (NCEP再解析のモデル格子データはこれ)。 さらに、高緯度では経度間隔を等しくすると距離間隔が狭くなるので、経度方向の点の数をまびくことがある (JRA55やJRA3Qのモデル格子はこれ)。 気象モデルに深入りしないデータ解析のためには、緯度・経度等間隔の格子に内挿してつかうことが多い。 (ただし緯度間隔と経度間隔はかならずしも同じではない。)

レーダー観測

オリジナルデータの位置は (距離、方位角、仰角) → 複数のレーダーをあわせてつかうときは緯度・経度・高さへ

衛星観測 (画像型センサー) のオリジナルデータ

オリジナルデータの位置は、衛星からみえた視野内の (x, y) → 地球の位置へは幾何変換が必要。

Registration: 格子点 (grid point) の値なのか、ます目 (grid box) の値なのか?

格子点データの計算機上でのもちかた

• 数字をふくむテキストファイルのこともある。
• 単純バイナリファイルのこともある。

数値の順序と地球上の位置との関係

• 数値の順序は、横書きの文書と同様に、東西方向は西から東へのことが多い。(出発点を0度にとる流儀、-180度にとる流儀がある。)
• 南北方向は、北から (横書き文書と同様)、南から (北緯が正の緯度がふえる方向) の両方がある。GrADS では南からを標準とし北からだと yrevと注記。
• 鉛直方向も、下から (zがふえる方向)、上から (pがふえる方向) の両方がある。

まわりかたは、(((東西), 南北), 上下) のことが多い。 Fortran の配列添え字は (lo, la, lev) となるが、Python では逆に [lev, la, lo] となる。

数値のバイナリ表現と圧縮

実数の値をとる数量の近似値は、32ビット (4バイト) または 64ビット (8バイト) の浮動小数点型であらわされることが多い。 仮数部の有効数字は、32ビット浮動小数点数ならば、21ビット程度。

実際の有効数字がもっとすくなければ、もっとすくないビット数で表現できるはず。

• ビット列としての圧縮 ... 浮動小数点型の格子データには、(0の連続を別にすれば) あまり効果がない。
• 数値としての圧縮
  • 同じ数値がつづくばあい 連長 (run length) 法
  • かぎられた種類の数値だけがあらわれるばあい、look up table 法 ... 表をつくって、各格子点については表内の番号をもつ。
  • scaling 法 ... 「実数値[格子点] = 整数値[格子点] × スケール値 + オフセット値」。各格子点については整数値をもつ。

データ欠損をあらわす値の問題

• 浮動小数点の IEEE 標準には NaN (not a number) をあらわすビットパタンがあるので、それがつかわれることも多い。
• 整数値にスケーリングしたばあいは整数値をつかった約束が必要。
• 欠損の理由を区別したいばあいもある。それぞれのデータセットごとのとりきめがある。

個別メタデータをどのようにもつか

メタデータ ... データについての情報。 ここでは臨時に、個別のファイルに何がふくまれているか (変数名、単位、日時、格子点番号と緯度経度の対応など) を「個別メタデータ」とよぶ。

単純バイナリファイルは、 個別メタデータ情報をファイル内にもたない。 かろうじてファイル名にふくめることができる。

個別メタデータ情報をファイル内にふくむ形式が提案されてきた。 それを解読するのには特殊なプログラムが必要になるが、よくつかわれる形式には対応できるプログラムが多く開発される。

GRIB (Grid Binary)

• 気象現業、とくに数値予報のための格子点データに特化した形式。WMOで制定された。
• version 1 と 2 がある。2 は GRIB2 ともいわれる。
  • (JRA55 は version 1、JRA3Q は GRIB2)。
  • 1 と 2 の考えかたは似ているが、直接の互換性はなく、別々のプログラムが必要である。 (Wesley Ebisuzaki さんの wgrib と wgrib2、GrADS の gribscan と grib2scan )。
• 複数のメタデータ部と格子点データ本体からなる。
• ビット単位のフィールドにデータをつめこむ。
• 気象要素名や単位などは、あらかじめ定義された番号による。 WMOに登録されたデータ提供センター (ECMWF, NCEP, JMA など) がそれぞれの定義表を提供している。
• 有効数字に必要最小限のビット数で数値を圧縮することができる。

わたしは、GrADSで単純バイナリに変換してそれをPython (numpy) でよむことが多い。

解析雨量と全国合成レーダーについては、ウェブ上にPythonで解読している例文をみつけたので、それを応用して自まえで numpy の配列に入れている。

NetCDF (CDFはもともとは common data format)

• 気象・海洋の研究者によくつかわれている形式。格子点データにかぎらない。
• UCAR (アメリカ合衆国の大気科学に関する大学の組合のようなもの) の Unidata プロジェクト (NCARとは横ならび) で制定された。
• NetCDF標準は、個別メタデータの形式まで。 個別メタデータの内容の意味は、CF 規約 (CF Convention) という約束にしたがうものが多い。 CF規約はインターネット上のコミュティの議論によって拡充されてきた。
• データ圧縮は、16ビット整数へのスケーリングはあるが、あまり柔軟ではない。

Python には netcdf4 というライブラリがつかわれてきたが、わたしは numpy との連携がしやすい xarray をつかっている。

[NetCDFとCF規約 (メモ 2010年)]

格子点データによる数量の2次元分布の表現

等値線

• matplotlib では contour。Cartopy 上でもつかえる。
  • オプションで、線の途中に値を数字ではさむことができる。

等値線間の段彩

• matplotlib では contourf。Cartopy 上でもつかえる。(GrADS の set gxout shaded に相当)
  • 色と数値との関係は colorbar をつかう。

格子点間隔がじゅうぶんこまかければ、等値線をひかず、画面を格子にくぎってそれぞれを色わけすればよい。 そのたぐいで、Cartopy で緯度経度で格子位置を指定してやる方法 (GrADS の set gxout grfill に相当) を、わたしはまだみつけておらず、 contourf にたよっている。 (imshow は画素の大きさが画面上で一定でないばあいには対応できないようだ。)

数値として内挿してはいけない定性的なデータの色分け表示

いまのところ、plot で marker をならべる方法

ベクトルの矢印表示

• matplotlib では quiver。Carotpy 上でもつかえる。ただしわたしはスケーリングをよく理解できていない。

風ベクトルの矢羽根表示

• matplotlib では barbs。Cartopy 上でもつかえる。

格子点でないデータの等値線

• (三角形分割して等値線のとおる位置をもとめる方法もあるが)
• 空間内挿 (scipy.interpolate.griddata linear) で格子点データをつくって contour にわたす

[空間内挿 (地理情報学教材 2000年)]

4次元時空間をいろいろな2次元で見る、フーリエ変換して振動数・波数領域で見る

• 廣田 勇, 1999: 気象解析学 -- 観測データの表現論。東京大学出版会, 175 pp. ISBN 4-13-060724-3. [読書メモ]

気象庁の地上観測の「記事」の雷の記録についてのコメント

定型文字列の処理。 (わたしが必要にせまられたら Awk でやるが、Python でもできるはず。 「自然言語処理」よりもビジネス応用が参考になるか?)

• 「.」を sep として split する
• 各文字列中で、数字4個の列をひろい (正規表現のマッチングでいけるはず)、時刻として解釈できるものの最初と最後の差をとる ... 継続時間と解釈できるだろう。 例外 (数字4個が時刻でないばあい、終了時刻がないばあいなど) のチェックは必要。
• ある日の記録が「－」でおわっているばあい、つぎの日の最初の「.」までをとりこむ。 (時刻は2400をたして解釈する。)