OCEAN COLOR: (NASA の衛星オープンデータ) の Python での取得と可視化
OCEAN COLOR WEB について
NASA は海洋観測のために、幾つかの衛星や装置でセンサーを保有しています。
これらのデータを統合して過去から現在までのアーカイブデータをオープンデータとして公開しています。
OCEAN COLOR とは
そのままではありますが、主に海洋での色のことです。
基本的には、海は太陽光の反射や散乱によってその色彩を変化させます。
しかしながら、それだけではなく海では様々な影響によってその色合い、色調を変えてゆきます。
例えば、プランクトンや含有金属、微生物などです。これらは自然の化学現象や人間の活動による結果だったりと地球規模での移り変わりを観測するための重要な指標の１つになっています。
海洋や気象にも密接に関係しているデータがこんなにも長期で獲得できるなんて素晴らしすぎますね。
今回は、身近な日本周辺のデータを取得して解析するために可視化したりして実際に動かし見ることで理解を深めるようなことをしていきたいと思います。
取得方法
今回は、 NASA の地球科学データ管理ツールの earthaccess を使用します。
この Python ライブラリでは、NASA EarthData のログイン情報を利用します。
以下で作成またはログインをして控えておきましょう。
Username:
Passward:
を使用します。
import earthaccess
auth = earthaccess.login(persist=True)
ここで先ほど控えた認証情報を入力します。
上記で一度認証すると、
!cat ~/.netrc
のように認証情報が保存されるのでこちらで手軽にするか、確認してみると良いでしょう。
これで OCEAN COLOR のデータセットの一覧を取得できます。
results = earthaccess.search_datasets(instrument="oci")
Datasets found: 20
このうちのレベル２の処理データを検索します。
results = earthaccess.search_data(
short_name="PACE_OCI_L2_AOP_NRT",
count=1,
)
Granules found: 4325
もう少し検索条件を絞って検索する場合は以下です。
tspan = ("2024-04-01", "2024-04-16")
bbox = (-76.75, 36.97, -75.74, 39.01)
clouds = (0, 50)
results = earthaccess.search_data(
short_name="PACE_OCI_L2_AOP_NRT",
temporal=tspan,
bounding_box=bbox,
cloud_cover=clouds,
)
Granules found: 2
では、１つ目のサンプルの情報を見てみます。
results[0]
上記のリンクを叩く事でも直接取得が可能です。
PACE_OCI.20240413T175656.L2.OC_AOP.V1_0_0.NRT.nc をダウンロードできます。
Python での取得する場合は以下です。
paths = earthaccess.open(results)
ここで取得するクラウドデータのパスを確認します。
paths
[<File-like object HTTPFileSystem, https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20240413T175656.L2.OC_AOP.V1_0_0.NRT.nc>,
<File-like object HTTPFileSystem, https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20240414T183158.L2.OC_AOP.V1_0_0.NRT.nc>]
次に、ローカルファイルにダウンロードする場合は以下です。
paths = earthaccess.download(results, "L2_AOP")
ダウンロードしたファイルの場所を確認します。
paths
['L2_AOP/PACE_OCI.20240413T175656.L2.OC_AOP.V1_0_0.NRT.nc',
'L2_AOP/PACE_OCI.20240414T183158.L2.OC_AOP.V1_0_0.NRT.nc']
ダウンロードするフォルダも指定することも可能です。
files = earthaccess.download(results, "./local_folder")
可視化
さてダウンロード方法もわかったので東京の衛星データを可視化してみましょう。
処理レベルが３つあるのでそれぞれの方法で記載していおります。
レベル１B
ここでは xarray という Python ライブラリを使用して処理レベル 1B の OCI L1B NetCDF ファイルを可視化します。
まずはエリアと時期を指定します。
tspan = ("2024-04-01", "2024-04-16")
bbox = (134.512,32.287, 141.576, 35.512)
この辺りを観測対象にしてます。
@ Google Map Satellite
その後、プロダクトを指定して検索します。
results = earthaccess.search_data(
short_name="PACE_OCI_L1B_SCI",
temporal=tspan,
bounding_box=bbox,
)
Granules found: 13
情報を取得してみます。
paths = earthaccess.open(results)
メタデータを確認します。
dataset = xr.open_dataset(paths[0])
dataset
ここでのメタデータは、観測値の一覧だけで実際に値を参照できないです。
なので実データ参照は以下です。
with h5netcdf.File(paths[0]) as file:
groups = list(file)
groups
['sensor_band_parameters',
'scan_line_attributes',
'geolocation_data',
'navigation_data',
'observation_data']
観測データにアクセスしてみます。
dataset = xr.open_dataset(paths[0], group="observation_data")
dataset
rhot_blue のデータの配列の形式を確認します。
dataset["rhot_blue"].sizes
Frozen({'blue_bands': 119, 'number_of_scans': 1710, 'ccd_pixels': 1272})
これでは形式が少しわかりにくいので図示します。
ブルーにも幾つかのバンドがあり、その中にスキャンとピクセルが縦と横の形式で格納されています。
バンド 10 を選択して可視化します。
plot = dataset["rhot_blue"].sel({"blue_bands": 10}).plot()
他にも red を選択して可視化もしてみます。
plot = dataset["rhot_red"].sel({"red_bands": 10}).plot()
レベル２
次に処理レベル２のプロダクトを処理します。
clouds = (0, 50)
results = earthaccess.search_data(
short_name="PACE_OCI_L2_AOP_NRT",
temporal=tspan,
bounding_box=bbox,
cloud_cover=clouds,
)
Granules found: 1
アクセスします。
paths = earthaccess.open(results)
レベル２では、geophysical_data が扱えるようになっています。
簡単にいうと、地球の形に沿って処理してある投影可能なデータです。
dataset = xr.open_dataset(paths[0], group="geophysical_data")
rrs = dataset["Rrs"]
rrs
では、可視化してみます。
plot = rrs.sel({"wavelength_3d": 100}).plot(cmap="viridis")
少し歪んでいるのが確認できると思います。日本の形もうっすら感じ取れそう、、、
この歪みを定量的なデータからよりわかりやすく可視化します。
そこで緯度経度を取得します。
dataset = xr.open_dataset(paths[0], group="navigation_data")
dataset = dataset.set_coords(("longitude", "latitude"))
dataset = dataset.rename({"pixel_control_points": "pixels_per_line"})
dataset = xr.merge((rrs, dataset.coords))
dataset
散布図で緯度経度の歪みを可視化します。
plot = dataset.sel(
{
"number_of_lines": slice(None, None, 1720 // 20),
"pixels_per_line": slice(None, None, 1272 // 20),
},
).plot.scatter(x="longitude", y="latitude")
先ほど得た位置座標を利用して再投影してみます。
rrs = dataset["Rrs"].sel({"wavelength_3d": 100})
plot = rrs.plot(x="longitude", y="latitude", cmap="viridis", vmin=0)
さらに、グリッドと海岸線を追加して地理空間データと融合して確認しやすくします。
fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines(draw_labels={"left": "y", "bottom": "x"})
plot = rrs.plot(x="longitude", y="latitude", cmap="viridis", vmin=0, ax=ax)
処理レベル２ならではの位置座標投影の情報があるからこその可視化や利用方法です。
位置座標的に局所的な変動を見る場合の分析もあります。
以下のように座標を限定します。
rrs_box = dataset["Rrs"].where(
(
(dataset["latitude"] > 36.52)
& (dataset["latitude"] < 36.55)
& (dataset["longitude"] < 142.46)
& (dataset["longitude"] > 142.43)
),
drop=True,
)
rrs_box.sizes
Frozen({'number_of_lines': 2, 'pixels_per_line': 3, 'wavelength_3d': 184})
限定された座標内での変化を可視化します。
rrs_stack = rrs_box.stack(
{"pixel": ["number_of_lines", "pixels_per_line"]}, create_index=False,
)
plot = rrs_stack.plot.line(hue="pixel")
レベル3
レベル 3 のファイル (L3M remote sensing reflectance (Rrs)) には、地球全体のグローバルマップが含まれています。
では、取得してみます。
tspan = ("2024-04-16", "2024-04-20")
results = earthaccess.search_data(
short_name="PACE_OCI_L3M_RRS_NRT",
temporal=tspan,
bounding_box=bbox,
)
Granules found: 30
データセットの確認をします。
dataset = xr.open_dataset(paths[0])
dataset
今回は、東南アジアからオセアニアにかけて可視化します。
fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines(draw_labels={"left": "y", "bottom": "x"})
rrs_442 = dataset["Rrs_442"].sel({"lat": slice(0, -60), "lon": slice(120, 180)})
plot = rrs_442.plot(cmap="viridis", vmin=0, ax=ax)
このように広域観測でも東京周辺を撮像していないこともあります。
なので期間を伸ばして、重ね合わせることでグローバルマップを生成します。
tspan = ("2024-04-12", "2024-04-24")
results = earthaccess.search_data(
short_name="PACE_OCI_L3M_CHL_NRT",
temporal=tspan,
granule_name="*.DAY.*.0p1deg.*",
)
Granules found: 13
一番近い日付から取得したデータセットを確認します。
dataset = xr.open_mfdataset(
paths,
combine="nested",
concat_dim="date",
)
dataset
取得したデータを可視化します。
chla = np.log10(dataset["chlor_a"])
chla.attrs.update(
{
"units": f'lg({dataset["chlor_a"].attrs["units"]})',
}
)
plot = chla.sel({"date": 0}).plot(aspect=2, size=4, cmap="GnBu_r")
さらに、雲の影響は NaN で抜け落ちているのでそれらを処理して平均値を可視化します。
chla_avg = chla.mean("date")
chla_avg.attrs.update(
{
"long_name": chla.attrs["long_name"],
"units": f'lg({chla.attrs["units"]})',
}
)
plot = chla_avg.plot(aspect=2, size=4, cmap="GnBu_r")
東京周辺も時系列的にも空間的にも補完ができて、日本列島が確認できると思います。
衛星比較
基本的には、光学衛星の観測になっています。
地球観測衛星の中でも広域の高頻度低解像度が多いです。
持続的に観測するために衛星の打ち上げ計画や現状の衛星も一覧であります。
Github
使用したコードは以下に保存しております。
環境構築
Python Version: 3.9 以上
pip install -r requirements.txt
earthacceess の環境構築は以下にもあります。
ノートブック
最後に
earthaccess を利用した NASA のオープンデータの取得と可視化を記載させて頂きました。
このようにデータ取得が簡潔になるような動きは多く、意味のある衛星データを統合して扱うからこそ見えてくるような解析が行えるような未来が見え始めてきて楽しいですね！！
おまけ
こちら以外にも記事執筆やコンペ解法記載をしているのでご参考になれば幸いです
衛星データ解析として、宙畑のライターもしています。
SAR 解析をよくやっていますが、画像系AI、地理空間や衛星データ、点群データに関心があります。
勉強している人は好きなので楽しく絡んでくれると嬉しいです。
Discussion