この記事の目次
ArgoデータをPythonへ読み込むだけなら数行で済みます。しかし、その数行だけでは、どのファイルを読むのか、生の値と補正済みの値のどちらを使うのか、どの品質管理フラグを受け入れるのか、鉛直座標が何を意味しているのかといった、適切なデータを選ぶための判断を意識しないまま、処理を進めてしまうことがあります。
これらを間違えても、たいていエラーにはならず、もっともらしいのに意図とは微妙に違う数値が出てくるだけです。そこでこの記事では、その選択肢を整理し、実際のフロートのデータでそれぞれの結果がどう変わるかを示したうえで、最初につまずきやすい箇所を挙げます。
Argoのファイルをまだ開いたことがない場合は、先に初心者向けArgo NetCDF形式の解説をご覧ください。この記事が前提とするArgo NetCDFファイルの構造を扱っています。
2つのやり方と、どこまで自分でやるか
ArgoのプロファイルをPythonで読む方法は、現実的には2つです。
argopyは、Euro-Argoの開発コミュニティが維持しているライブラリで、検索、ダウンロード、キャッシュ、品質による絞り込みまでを引き受け、xarrayのDatasetを返します。xarrayとnetCDF4の組み合わせは、GDACのファイルを直接読み、その中身をそのまま渡します。
| argopy | xarray＋生NetCDF | |
|---|---|---|
| ファイルの特定 | 海域・フロート・プロファイルで検索 | 自分でパスを特定する |
| QCとデータモード | ユーザーモードに応じて自動適用 | 自分で適用する |
| 出力の形 | 整えられたDataset、再整形も可能 | ファイル本来の次元 |
| 絞り込み前のデータ | expertモードで取得できる | 常にそのまま得られる |
生のNetCDFを読むほうが上級者向け、というわけではありません。ある海域や特定のフロートを調べたい、というところから始まる解析であれば、まずargopyを使うのが無難です。生のNetCDFが適しているのは、argopyが公開していないものが必要なとき、ローカルミラーで作業しているとき、ファイルに何の加工も入っていない状態の中身を正確に知りたいときです。
argopyなら3行
from argopy import DataFetcher
ds = DataFetcher(src="erddap", mode="standard").float(6902746).load().data返ってくるのは「点」の形をしたxarrayのDatasetです。すべてのサイクルのすべての測定点がN_POINTS次元に沿った1行になり、LATITUDE、LONGITUDE、TIMEが座標として付きます。絞り込みには便利ですがプロファイルの描画には向かないので、argopyにはこれを並べ替えるメソッドがあります。
profiles = ds.argo.point2profile()これで見慣れたN_PROF × N_LEVELSの形になります。
最初の1行のうち、2つの引数がほかより重要です。
srcはデータソースを選びます。erddap（デフォルト）、gdac（公式サーバー、またはGDAC互換のローカルディレクトリ）、argovisのいずれかです。生物地球化学変数を取得できるのはerddapとgdacで、argovisはコアデータのみです。modeは、argopyがどこまで処理した状態のデータを返すかを選びます。
ユーザーモードの選択は解析結果を左右する
argopyには3つのユーザーモードがあり、どれを選ぶかでDatasetに含まれる測定値そのものが変わります。
expertモードは、データソースが保持しているままのデータを絞り込みなしで返します。生の変数と補正済み変数の両方があり、すべての品質管理フラグが含まれます。standardモード（既定）は、品質管理フラグが1または2の測定値を残し、PARAM_ADJUSTEDをPARAMへ統合して補正済み変数を削除し、専門家向けの変数を隠します。researchモードはさらに厳しく、遅延モードのデータのみ、品質管理フラグは1のみ、加えてコア変数には圧力誤差の上限が課されます。
138サイクルの遅延モードのフロート、WMO 6902746について、3つのモードが返すものは次のとおりです。
| モード | 測定値 | プロファイル | 変数 |
|---|---|---|---|
expert | 14,345 | 138 | 23 |
standard | 12,518 | 118 | 15 |
research | 12,518 | 118 | 9 |
expertモードの結果は、GDACにあるこのフロートのマルチプロファイルファイルと完全に一致します。取得元のデータがそのまま渡されていることは、ここで確かめられます。一方standardモードでは、測定値の約13%とプロファイル20本が落ちています。argopyが誤って捨てているわけではなく、Argo計画がすでに下した品質判断がそのまま適用されているだけです。
スクリプトが出したプロファイル数は、そのフロートの情報であると同時に、自分がどう絞り込んだかの情報でもあります。モードは結果とあわせて記録してください。
argopyが適用する前のフラグとデータモードの意味は、Argoデータの品質管理が参考になります。
生NetCDFを直接読む
GDACのファイルを直接読む場合も、コード自体は1行です。
import xarray as xr
ds = xr.open_dataset("6902746_prof.nc")返ってくるのは、ファイルがもともと持っている形そのままです。このフロートではN_PROFが138、N_LEVELSが110で、ほかにN_PARAM、N_CALIB、N_HISTORYがあります。絞り込みは何もかかっておらず、TEMPとTEMP_ADJUSTEDが両方入っているので、どちらを使うかは自分で決めることになります。
測定値の配列は長方形ですが、実際に入っているプロファイルの長さは揃っていません。すべてのプロファイルがファイル中でもっとも長いものに合わせて埋められるため、短いプロファイルのあとには意味を持たない欠損値が続きます。
このファイルでは、ほとんどのプロファイルが100を超える測定点を埋めていますが、なかにはずっと少ないものもあり、最も短いプロファイルは16点しかありません。結果として、配列全体の約5%が中身のないパディングになっています。N_LEVELSの方向に平均などの統計量を計算するとき、この部分を除外しないと、中身のない値まで一緒に集計してしまいます。
最初に間違えやすい4つのこと
はじめの3つが問題になるのは、生のNetCDFを直接読む場合です。argopyはフラグとデータモードの型変換を済ませたうえで、expertモード以外では生の値と補正済み変数の選択も代わりに行います。最後の1つは、どちらの方法で読んでも同じように効いてきます。
品質管理フラグは整数ではなくバイト列
QCの変数は1文字として保存されているため、xarrayが返すのはバイト列を要素に持つobject配列で、値がない場所には浮動小数点のNaNが入ります。この配列を整数の1と比較しても一致せず、並べ替えようとするとバイト列と浮動小数点を比較できないためTypeErrorになります。バイトリテラルと比較してください。
import numpy as np
good = np.isin(ds["TEMP_QC"].values, [b"1", b"2"])DATA_MODEにも同じ問題があります。値は想定しがちな文字列ではなく、b"R"、b"A"、b"D"です。
圧力は深さではない
Argoの鉛直座標はdbar単位の圧力であり、mとの数値の近さは目安にすぎません。重力は緯度によって変わるため、正しく換算するには緯度が必要です。
import gsw
depth_m = -gsw.z_from_p(2000.0, 30.0)2,000 dbarでは、緯度30°で約1,977 m、赤道でおよそ1,980 m、緯度60°で約1,972 mになります。差は1%未満ですが、ランダムな誤差ではなく、同じ緯度なら必ず同じだけずれます。そのため、データを重ねても打ち消されません。しかもずれの大きさは緯度によって変わるため、離れた海域のプロファイルを比べると、その差が見かけ上の違いとして残ります。しきい値の計算であれば、dbarのまま扱えばこの問題自体が起きません。
生の値と補正済み変数を混ぜない
TEMPとTEMP_ADJUSTED_QCを組み合わせたり、補正済みの水温と生の塩分を並べたりすると、どちらの状態にも対応しないプロファイルができてしまいます。まずデータモードを読み、対応する値とフラグの組を選び、すべての変数で組み合わせを揃えてください。生物地球化学変数では、モードは変数ごとです。DATA_MODEではなくPARAMETER_DATA_MODEを見ます。
単位は変数の属性に書かれていて、想定どおりとは限らない
PRESにはdecibar、TEMPにはdegree_Celsius、PSALにはpsuが付いています。psuは実用塩分であり、絶対塩分ではありません。密度や保存水温を導くには、まずTEOS-10による変換が必要です。溶存酸素はDOXYとしてµmol/kgで届き、センサーの出力からすでに換算されています。別のデータセットと組み合わせるときにこれが効いてくる理由は、海洋データの溶存酸素の単位で扱っています。
解析環境は先に固定しておく
Argo関連のツールは、科学計算向けPythonの多数のライブラリに依存しており、動作するバージョンの組み合わせが必ずしも「すべて最新」とは限りません。そのため、今日は動く解析が半年後に作り直した環境では動かないこともありますが、環境の記録を残しておけば避けられます。
正確なバージョンをlockfile、環境定義ファイル、あるいはコンテナイメージとして記録しておけば、「自分の環境では動く」という状態を、査読者や将来の自分が再現できる状態に変えられます。手間をかけずにこれを行う方法は研究の解析環境を再現可能にする方法：Docker入門で扱っています。
先に見ておくとコードは短くなる
Argoを扱い始めた最初の1週間に書くコードの多くは、解析そのものではなく、データがどうなっているかを知るためのものです。具体的には、その海域にどのフロートがいるのか、あるフロートのサイクルは対象期間をカバーしているのか、酸素は入っているのか、探している特徴がそもそもあるのか、といった点を確認します。
これはデータを読む作業なので、読み込み処理を書いて確かめるより、目で見たほうが速く済みます。どのWMO IDとサイクルが重要か分かれば、argopyの.float()や.profile()で対象のフロートとプロファイルを直接取得でき、探す処理を省いて解析のコードに集中できます。
こうしてOceanGraph上でデータを事前に絞り込む作業については、Pythonを使わずにArgoフロートデータを可視化する手順で扱っており、App GuideにはSearch and Bookmark（英語）があります。いずれもPythonを不要にするものではなく、何を調べたいかがはっきりしないうちからPythonを書き始めずに済むようにするものです。
よくある質問
argopyと生のNetCDFのどちらを使うべきですか
理由がないかぎりargopyを使ってください。検索、キャッシュ、品質の判断を一貫して処理してくれますし、絞り込まれていないデータが必要なときはexpertモードがあります。ファイルを直接読むのは、ローカルミラーで作業しているとき、argopyが公開していない変数が必要なとき、ファイルの内容を正確に確認したいときです。
モードを変えるとプロファイル数が変わるのはなぜですか
モードごとに適用される品質管理とデータモードの絞り込みが違うためです。standardモードはフラグ1と2を残し、researchモードは遅延モードでフラグ1のものだけを残します。絞り込みの結果、残る測定値が一つもなくなったプロファイルは、Datasetに含まれなくなります。
品質管理フラグで自分で絞り込むには
どちらの方法で読んだかで変わります。ファイルを直接読む場合は、フラグが数値ではなく文字として保存されているので、b"1"やb"2"のようなバイトリテラルと比較します。argopyはQCの変数を整数へ変換して返すため、expertモードのDatasetを絞り込むときはds.argo.filter_qc(QC_list=[1, 2])のように整数で指定します。
gswなしで圧力を深さへ換算できますか
近似はできますが、緯度を考慮した換算はgswの関数1つで済み、避けられる誤差を取り除けます。解析にmが必要でないなら、dbarのまま扱うほうが簡単ですし、Argoがもともとdbarで配布していることとも合っています。
最初に見るだけならPythonは必要ですか
必要ありません。むしろ、コードを書く前にプロファイルを目で確かめておくほうが速く、コードで処理すべき内容も絞れることがあります。Pythonが必要になるのは、繰り返し実行できる処理や独自の導出量が要るとき、あるいは1つずつ開いていられないほど多くのプロファイルを扱うときです。
