NASA Earthdata の衛星データを使って Livebook で CO2(二酸化炭素)の大気濃度を可視化する

はじめに

SDGs の 13 番目の目標に「気候変動に具体的な対策を」というのがあります

近年では二酸化炭素などの温室効果ガスを削減するよう、あちこちで言われていますね

「人間が気候変動を引き起こしていることはどのようにしてわかるか？」については @zacky1972 さんの記事を参考してください

では、具体的にどれくらい大気中の二酸化炭素濃度は上昇しているのでしょうか

というわけで、今回は NASA が出している OCO-2 という衛星のデータを使って、大気中の二酸化炭素濃度を可視化します

OCO-2 の公式サイトはこちら

全て英語ですが、 OCO-2 がどのようにして二酸化炭素濃度を取得しているのかなど、動画も使って解説しています

また、今回は日本の Tellus ではなく、 NASA の Earthdata からデータをダウンロードします

Earthdata

今回も可視化には Elixir の Livebook を使っていきます

今回実装したノートブックはこちら

事前準備

Earthdata は全て英語ですが、かなりの量のデータが公開されています

これを使わない手はありません

データのダウンロードにはユーザー登録が必要になります

以下のサイトに日本語で Earthdata のユーザー登録の手順を説明してくれています（私もここを見て登録しました）

データの取得

Earthdata は様々な項目で衛星データを検索できるようになっています

以下のリンクは OCO-2 公式サイトにも貼られている、 OCO-2 データの検索結果です

OCO-2 のデータにも色々ありますが、今回は　OCO2_L2_Lite_FP 11r （以下の画像内の先頭に表示されているもの）を使います

衛星データは「どれくらい加工されたか」でレベル0から4までに分類されており、今回使うのはレベル２のデータです

衛星のセンサーが読み取った値そのものではなく、それを目的のために扱いやすくしたものです

今回でいうと二酸化炭素濃度が値としてそのまま取得できます

また、 OCO2_L2_Lite では一定の精度以上で取得できた値だけに絞り込まれています

OCO2_L2_Lite_FP 11r の Subset / Get Data をクリックしましょう

そうすると以下のようなモーダルが開いたと思います

ここで場所や時間を指定してデータを絞り込めるようになっています

OCO2_L2_Lite_FP 11r は1ファイルに地球全体の1日分のデータが入っているため、場所の絞り込みは不要です

今回は12月中のデータを見たいので、何度か検索して以下のデータをダウンロードします

• 2016年12月1日
• 2017年12月1日
• 2018年12月1日
• 2019年12月1日
• 2020年12月1日
• 2021年12月1日
• 2022年12月1日 - 2022年12月31日 (1ヶ月分)

ダウンロードしたファイルは以下のようなファイル名になっています

• oco2_LtCO2_161201_B11014Ar_230111183332s.nc4
• oco2_LtCO2_221201_B11014Ar_230118175544s.nc4

これを以下のように年月毎のディレクトリーに格納したものとして進めます

• tmp
  • OCO2
    • 201612
      • oco2_LtCO2_161201_B11014Ar_230111183332s.nc4
    • 201712
      • oco2_LtCO2_171201_B11014Ar_221129005324s.nc4
    • 201812
      • oco2_LtCO2_181201_B11014Ar_221017172333s.nc4
    • ...
    • 202212
      • oco2_LtCO2_221201_B11014Ar_230118175544s.nc4
      • oco2_LtCO2_221202_B11014Ar_230118180140s.nc4
      • oco2_LtCO2_221203_B11014Ar_230118180210s.nc4
      • ...

実行環境

• Elixir: 1.14.2 OTP 24
• Livebook: 0.8.1

以下のリポジトリーの Docker コンテナ上で起動しています

Docker が使える環境であれば簡単に実行できます

https://docs.docker.com/engine/install/

Docker Desktop を無償利用できない場合は Rancher Desktop を使ってください

https://rancherdesktop.io/

今回 NetCDF という、気象などのデータを扱うモジュールをインストールするため、 Dockerfile に Rust のインストールを追加しています

Elixir でデータ解析をしようという人が増えている証拠ですね

https://github.com/DockYard/netcdf

Elixir 版 NetCDF は Rust 版 NetCDF を NIF で呼び出しています

そのため、ビルド用に Rust をインストールする必要がありました

https://github.com/georust/netcdf

すでに私の elixir-learning リポジトリーをクローンしている人は Pull してください

私のリポジトリーを使う場合、以下の手順で Livebook を起動できます

git clone https://github.com/RyoWakabayashi/elixir-learning.git

cd elixir-learning

$ docker-compose up
...
Attaching to livebook
livebook    | [Livebook] Application running at http://0.0.0.0:<ポート番号>/?token=<認証トークン>

最後に表示される URL をブラウザで開けば Livebook にアクセスできます

右上 New notebook をクリックし、新しいノートブックを開けば準備完了です

セットアップ

Mix.install(
  [
    {:nx, "~> 0.4"},
    {:evision, "~> 0.1"},
    {:exla, "~> 0.4"},
    {:netcdf, "~> 0.1"},
    {:geo, "~> 3.4"},
    {:kino, "~> 0.8"},
    {:kino_maplibre, "~> 0.1"},
    {:kino_vega_lite, "~> 0.1"}
  ],
  config: [
    nx: [
      default_backend: EXLA.Backend
    ]
  ]
)

• Nx: 行列演算
• Evision: 画像処理
• EXLA: 行列演算の高速化
• NetCDF: NetCDF 形式データの入出力
• Geo: GeoJSON 形式データの入出力
• Kino: Livebook での入出力
• Kino MapLibre: 地図出力
• Kino VegaLite: グラフ出力

また、 Nx による行列演算を EXLA を使って高速化する設定をしています

衛星データの読込

衛星データを格納したディレクトリーを指定します

oco2_dir = "/tmp/OCO2/"

まず最初に2022年12月のデータを対象として処理していきます

target_month_dir = "202212"

対象ディレクトリー内の最初のファイル（2022年12月1日分）を開きます

{:ok, file} =
  [oco2_dir, target_month_dir]
  |> Path.join()
  |> File.ls!()
  |> Enum.sort()
  |> Enum.at(0)
  |> then(&[oco2_dir, target_month_dir, &1])
  |> Path.join()
  |> NetCDF.File.open()

実行結果は以下のようになります

{:ok,
 %NetCDF.File{
   resource: #Reference<0.2046660988.1870790657.260263>,
   filename: "/tmp/OCO2/202212/oco2_LtCO2_221201_B11014Ar_230118175544s.nc4",
   variables: ["sounding_id", "levels", "bands", "vertices", "footprints", "date", "latitude",
    "longitude", "time", "solar_zenith_angle", "sensor_zenith_angle", "xco2_quality_flag",
    "xco2_qf_bitflag", "xco2_qf_simple_bitflag", "source_files", "file_index", "vertex_latitude",
    "vertex_longitude", "xco2", "xco2_uncertainty", "xco2_apriori", "pressure_levels",
    "co2_profile_apriori", "xco2_averaging_kernel", "pressure_weight"]
 }}

variables に NetCDF が持っているデータの項目名が入っています

NetCDF.Variable.load で項目毎にデータを読み込むことができます

今回は以下の５つの項目を使います

• latitude (緯度)
• longitude　(経度)
• vertex_latitude (頂点の緯度)
• vertex_longitude (頂点の経度)
• xco2 (二酸化炭素濃度)

variables =
  ["latitude", "longitude", "vertex_latitude", "vertex_longitude", "xco2"]
  |> Enum.map(fn var_name ->
    file
    |> NetCDF.Variable.load(var_name)
    |> elem(1)
  end)

実行結果は以下のようになります

[
  %NetCDF.Variable{
    name: "latitude",
    value: [-43.986419677734375, -43.99115753173828, -43.99583053588867, -44.00041961669922,
    ...
     -43.913028717041016, -43.917449951171875, -43.92179489135742, ...],
    type: :f32,
    attributes: %{
      "comment" => ["center latitude of the measurement"],
      "long_name" => ["latitude"],
      "missing_value" => -999999.0,
      "units" => ["degrees_north"]
    }
  },
  %NetCDF.Variable{
    name: "longitude",
    value: [-151.19435119628906, -151.21084594726562, -151.22750854492188, -151.24432373046875,
    ...
      "missing_value" => -999999.0,
      "units" => ["ppm"]
    }
  }
]

項目毎に名前、値（配列）、型、単位などの情報が格納されています

配列のままだと計算しにくいことがあるので、テンソルに変換します

as_nx_type = fn
  :i8 -> :s8
  :i16 -> :s16
  :i32 -> :s32
  :i64 -> :s64
  t -> t
end

var_tensors =
  for var <- variables, into: %{} do
    tensor = Nx.tensor(var.value, type: as_nx_type.(var.type))
    {var.name, tensor}
  end

実行結果は以下のようになります

%{
  "latitude" => #Nx.Tensor<
    f32[154515]
    EXLA.Backend<host:0, 0.921923193.662831120.257644>
    [-43.986419677734375, -43.99115753173828, -43.99583053588867, -44.00041961669922, ...]
  >,
  "longitude" => #Nx.Tensor<
    f32[154515]
    EXLA.Backend<host:0, 0.921923193.662831120.257645>
    [-151.19435119628906, -151.21084594726562, -151.22750854492188, -151.24432373046875, ...]
  >,
  "vertex_latitude" => #Nx.Tensor<
    f32[618060]
    EXLA.Backend<host:0, 0.921923193.662831120.257646>
    [-43.9932975769043, -43.974761962890625, -43.979583740234375, -43.99810791015625, ...]
  >,
  "vertex_longitude" => #Nx.Tensor<
    f32[618060]
    EXLA.Backend<host:0, 0.921923193.662831120.257647>
    [-151.18304443359375, -151.18838500976562, -151.2056427001953, -151.20030212402344, ...]
  >,
  "xco2" => #Nx.Tensor<
    f32[154515]
    EXLA.Backend<host:0, 0.921923193.662831120.257648>
    [414.4909973144531, 415.3614807128906, 414.925048828125, 414.541015625, 414.994140625, ...]
  >
}

衛星データのグラフ表示

ざっくりデータの傾向を見てみるためにグラフ化してみます

display_graph = fn tensors, item_name ->
  plot_data =
    tensors[item_name]
    |> Nx.to_flat_list()
    |> Enum.with_index()
    |> Enum.map(fn {item, index} ->
      %{
        "index" => index,
        item_name => item
      }
    end)

  VegaLite.new(width: 700)
  |> VegaLite.data_from_values(plot_data)
  |> VegaLite.mark(:line)
  |> VegaLite.encode_field(:x, "index", type: :quantitative)
  |> VegaLite.encode_field(:y, item_name, type: :quantitative)
end

まず経度をグラフ化してみます

display_graph.(var_tensors, "longitude")

-180 から 180 の間を1日かけて動いているのがわかります（時々ヒゲのように出ているところはありますが）

これは地球の自転によって衛星の観測場所が移動していることを表しています

続いて緯度をグラフ化してみます

display_graph.(var_tensors, "latitude")

-90 から 60 あたりを1日に何度も行ったり来たりしています

これは衛星が南北方向にグルグル回っている様子です

最後に二酸化炭素濃度をグラフ化します

display_graph.(var_tensors, "xco2")

400 の少し上あたりで小刻みに動いています

これだと線が上の方に集中していてみにくいので、表示範囲を最小値と最大値の間に限定します

{max_co2, min_co2} =
  {
    var_tensors["xco2"] |> Nx.reduce_max() |> Nx.to_number(),
    var_tensors["xco2"] |> Nx.reduce_min() |> Nx.to_number(),
  }

plot_data =
  var_tensors["xco2"]
  |> Nx.to_flat_list()
  |> Enum.with_index()
  |> Enum.map(fn {item, index} ->
    %{
      "index" => index,
      "xco2" => item
    }
  end)

y_scale = %{"domain" => [min_co2, max_co2]}

VegaLite.new(width: 700)
|> VegaLite.data_from_values(plot_data)
|> VegaLite.mark(:line)
|> VegaLite.encode_field(:x, "index", type: :quantitative)
|> VegaLite.encode_field(:y, "xco2", type: :quantitative, scale: y_scale)

二酸化炭素濃度が地球の全域でおよそ 418 ppm くらいであることが分かります

この緯度、経度、二酸化炭素濃度のデータ数は同じになっていて、これらを地図上にプロットすれば、地球上のどこがどの程度の二酸化炭素濃度なのか可視化できそうです

衛星データの地図表示

地図に表示するにあたって、二酸化炭素濃度を色として表したいと思います

色は 0 から 255 の範囲で表現するため、現在の値をその範囲内に収まるよう補正します

corrected_co2 =
  var_tensors["xco2"]
  |> Nx.subtract(min_co2)
  |> Nx.multiply(255 / (max_co2 - min_co2))
  |> Nx.as_type(:u8)

黒から白の色変化では分かりづらいので、ヒートマップを準備します

# Evision(OpenCV)のJETヒートマップを準備
jet_map =
  Nx.iota({256})
  |> Nx.tile([256, 1])
  |> Nx.as_type(:u8)
  |> Evision.Mat.from_nx_2d()
  |> Evision.applyColorMap(Evision.Constant.cv_COLORMAP_JET())

JETヒートマップでは値が小さい方から順に青->水->緑->黄->赤に変化します

二酸化炭素が多いほど暑いイメージと合致するため、このヒートマップを使えば直感的に伝わりそうです

MapLibre でヒートマップ表現するためのオプションを生成します

common_expression =
  ["step", ["get", "xco2"]]

bgr =
  jet_map
  |> Evision.Mat.to_nx(EXLA.Backend)
  |> then(& &1[[0, 0..255]])

b = bgr[[0..255, 0]]
g = bgr[[0..255, 1]]
r = bgr[[0..255, 2]]

get_expression = fn table_tensor ->
  [_ | tail] =
    [Nx.iota({256}), table_tensor]
    |> Nx.stack()
    |> Nx.transpose()
    |> Nx.to_flat_list()
  
    tail
end

r_expression = common_expression ++ get_expression.(r)
g_expression = common_expression ++ get_expression.(g)
b_expression = common_expression ++ get_expression.(b)

jet_expression =
  [
    "rgb",
    ["round", r_expression],
    ["round", g_expression],
    ["round", b_expression],
  ]

衛星データを MapLibre に読み込ませるための形式にします

map_source =
  %{
    "latitude" => Nx.to_flat_list(var_tensors["latitude"]),
    "longitude" => Nx.to_flat_list(var_tensors["longitude"]),
    "xco2" => Nx.to_flat_list(corrected_co2)
  }

地図に表示します

データがあるところに二酸化炭素濃度に応じた色の丸を置いていきます

MapLibre.new()
|> MapLibre.add_table_source(
  "xco2",
  map_source,
  {:lng_lat, ["longitude", "latitude"]},
  properties: ["xco2"]
)
|> MapLibre.add_layer(
  id: "xco2_layer",
  source: "xco2",
  type: :circle,
  paint: [circle_color: jet_expression]
)

チンアナゴのようにニョロニョロとした縦の線が何本も見えます

これが OCO-2 の衛星軌道です

何度も南北に縦断しながら東西に移動しているのが分かると思います

南極の上空に、二酸化炭素濃度が非常に低い場所、非常に高い場所がそれぞれあり、大体の場所は黄緑からオレンジ（中間の値）になっています

また、観測できるのはあくまでも衛星が通ったところだけなので、全ての座標で二酸化炭素濃度が得られるわけではなさそうです

精度が低かったところはデータから省かれているため、そういった場所の二酸化炭素濃度も得られません

1ヶ月分の衛星データを表示する

1日分だと足りないかもしれないので、1ヶ月分を見てみます

ただし地球全域にすると多すぎるので、日本近辺に限定します

load_japan_data = fn filename ->
  {:ok, nc_file} =
    [oco2_dir, target_month_dir, filename]
    |> Path.join()
    |> NetCDF.File.open()

  longitude =
    nc_file
    |> NetCDF.Variable.load("longitude")
    |> elem(1)
    |> then(&Nx.tensor(&1.value, type: as_nx_type.(&1.type)))

  latitude =
    nc_file
    |> NetCDF.Variable.load("latitude")
    |> elem(1)
    |> then(&Nx.tensor(&1.value, type: as_nx_type.(&1.type)))

  is_japan =
    # 経度127-147,緯度26-46の範囲を対象にする
    Nx.logical_and(
      Nx.logical_and(Nx.greater_equal(longitude, 127), Nx.less(longitude, 147)),
      Nx.logical_and(Nx.greater_equal(latitude, 26), Nx.less(latitude, 46))
    )

  japan_indeces =
    is_japan
    |> Nx.to_flat_list()
    |> Enum.with_index()
    |> Enum.filter(fn {is_true, _} -> is_true > 0 end)
    |> Enum.map(fn {_, index} -> index end)

  xco2 =
    nc_file
    |> NetCDF.Variable.load("xco2")
    |> elem(1)
    |> then(&Nx.tensor(&1.value, type: as_nx_type.(&1.type)))

  case japan_indeces do
    [] ->
      nil
    japan_indeces ->
      [longitude, latitude, xco2]
      |> Nx.stack()
      |> Nx.take(Nx.tensor(japan_indeces), axis: 1)
  end
end

対象月のディレクトリーにある全ファイルを読み込み、結果を結合します

target_month_tensor =
  [oco2_dir, target_month_dir]
  |> Path.join()
  |> File.ls!()
  |> Enum.sort()
  |> Enum.map(&load_japan_data.(&1))
  |> Enum.filter(& &1 != nil)
  |> Nx.concatenate(axis: 1)

対象範囲内の二酸化炭素濃度の最小値、最大値を算出し、色を補正します

{max_co2, min_co2} =
  {
    target_month_tensor[2] |> Nx.reduce_max() |> Nx.to_number(),
    target_month_tensor[2] |> Nx.reduce_min() |> Nx.to_number(),
  }

target_month_corrected_co2 =
  target_month_tensor[2]
  |> Nx.subtract(min_co2)
  |> Nx.multiply(255 / (max_co2 - min_co2))
  |> Nx.as_type(:u8)

map_source =
  %{
    "latitude" => Nx.to_flat_list(target_month_tensor[1]),
    "longitude" => Nx.to_flat_list(target_month_tensor[0]),
    "xco2" => Nx.to_flat_list(target_month_corrected_co2)
  }

MapLibre.new(center: {136, 36}, zoom: 3)
|> MapLibre.add_table_source(
  "xco2",
  map_source,
  {:lng_lat, ["longitude", "latitude"]},
  properties: ["xco2"]
)
|> MapLibre.add_layer(
  id: "xco2_layer",
  source: "xco2",
  type: :circle,
  paint: [circle_color: jet_expression]
)

1ヶ月分を使っても、日本の各都道府県単位でカバーすることはできなさそうです

多角形による地図表示

OCO-2 のデータには vertex_latitude や vertex_longitude がありました

これは二酸化炭素濃度の観測範囲を４つの頂点で表したものです

なので latitude の配列の長さ * 4 = vertex_latitude の配列の長さになっています

先ほどまでは適当に丸を付けていましたが、実際にはこの領域を塗りつぶすのが正しいようです

全ての観測領域を個別に計算すると時間がかかるので、また日本近辺だけにします

polygon_tensor =
  [
    var_tensors["vertex_longitude"],
    var_tensors["vertex_latitude"]
  ]
  |> Nx.stack()
  |> Nx.transpose()
  # 観測範囲数 * 4(頂点数) * 2(経度と緯度)にする
  |> Nx.reshape({:auto, 4, 2})

co2_polygon_tensor =
  [
    polygon_tensor,
    corrected_co2 |> Nx.new_axis(-1) |> Nx.tile([4]) |> Nx.new_axis(-1)
  ]
  |> Nx.concatenate(axis: 2)

is_japan =
  Nx.logical_and(
    Nx.logical_and(
      Nx.greater_equal(var_tensors["longitude"], 127), 
      Nx.less(var_tensors["longitude"], 147)
    ),
    Nx.logical_and(
      Nx.greater_equal(var_tensors["latitude"], 26), 
      Nx.less(var_tensors["latitude"], 46)
    )
  )

japan_indeces =
  is_japan
  |> Nx.to_flat_list()
  |> Enum.with_index()
  |> Enum.filter(fn {is_true, _} -> is_true > 0 end)
  |> Enum.map(fn {_, index} -> index end)
  |> Nx.tensor()

japan_co2_polygon_tensor = Nx.take(co2_polygon_tensor, japan_indeces, axis: 0)

この時点の実行結果は以下のようになっています

#Nx.Tensor<
  f32[366][4][3]
  EXLA.Backend<host:0, 0.921923193.662831115.258800>
  [
    [
      [142.25299072265625, 26.05797576904297, 179.0],
      [142.24713134765625, 26.075912475585938, 179.0],
      [142.25540161132812, 26.063196182250977, 179.0],
      [142.26123046875, 26.045265197753906, 179.0]
    ],
    [
      [142.25537109375, 26.06390953063965, 178.0],

366(1日分の日本近辺の観測領域数) * 4(頂点数) * 3(経度、緯度、色) のテンソルになりました

この多角形の配列を GeoJSON に変換します

geo_collection =
  japan_co2_polygon_tensor
  |> Nx.to_batched(1)
  |> Enum.map(fn co2_polygon ->
    coordinates =
      co2_polygon[[0, 0..3, 0..1]]
      |> Nx.to_flat_list()
      |> Enum.chunk_every(2)
      |> Enum.map(&List.to_tuple(&1))

    co2 =
      co2_polygon[[0, 0, 2]]
      |> Nx.to_number()

    %Geo.Polygon{
      coordinates: [coordinates],
      properties: %{"xco2" => co2}
    }
  end)
  |> Enum.to_list()
  |> then(& %Geo.GeometryCollection{geometries: &1})

実行結果は以下のようになります

多角形毎に一つの色（二酸化炭素濃度の相対値）が入っています

%Geo.GeometryCollection{
  geometries: [
    %Geo.Polygon{
      coordinates: [
        [
          {142.25299072265625, 26.05797576904297},
          {142.24713134765625, 26.075912475585938},
          {142.25540161132812, 26.063196182250977},
          {142.26123046875, 26.045265197753906}
        ]
      ],
      srid: nil,
      properties: %{"xco2" => 179.0}
    },
    %Geo.Polygon{
...
}

JETカラーマップで地図にプロットします

MapLibre.new(center: {141, 30}, zoom: 4)
|> MapLibre.add_geo_source("co2", geo_collection)
|> MapLibre.add_layer(id: "co2_polygon", source: "co2", type: :fill, paint: [fill_color: jet_expression])

すごく小さいですが、ズームしていくと平行四辺形が並んでいるのが分かります

観測領域は地球規模で見るとかなり小さく、これで地球全域をカバーするのは無理だな、ということが分かります

でもせっかくなので、これも1ヶ月分表示してみます

load_japan_polygon_data = fn filename ->
  {:ok, nc_file} =
    [oco2_dir, target_month_dir, filename]
    |> Path.join()
    |> NetCDF.File.open()

  longitude =
    nc_file
    |> NetCDF.Variable.load("longitude")
    |> elem(1)
    |> then(&Nx.tensor(&1.value, type: as_nx_type.(&1.type)))

  latitude =
    nc_file
    |> NetCDF.Variable.load("latitude")
    |> elem(1)
    |> then(&Nx.tensor(&1.value, type: as_nx_type.(&1.type)))

  is_japan =
    Nx.logical_and(
      Nx.logical_and(Nx.greater_equal(longitude, 127), Nx.less(longitude, 147)),
      Nx.logical_and(Nx.greater_equal(latitude, 26), Nx.less(latitude, 46))
    )

  japan_indeces =
    is_japan
    |> Nx.to_flat_list()
    |> Enum.with_index()
    |> Enum.filter(fn {is_true, _} -> is_true > 0 end)
    |> Enum.map(fn {_, index} -> index end)

  xco2 =
    nc_file
    |> NetCDF.Variable.load("xco2")
    |> elem(1)
    |> then(&Nx.tensor(&1.value, type: as_nx_type.(&1.type)))

  vertex_longitude =
    nc_file
    |> NetCDF.Variable.load("vertex_longitude")
    |> elem(1)
    |> then(&Nx.tensor(&1.value, type: as_nx_type.(&1.type)))

  vertex_latitude =
    nc_file
    |> NetCDF.Variable.load("vertex_latitude")
    |> elem(1)
    |> then(&Nx.tensor(&1.value, type: as_nx_type.(&1.type)))

  polygon_tensor =
    [
      vertex_longitude,
      vertex_latitude
    ]
    |> Nx.stack()
    |> Nx.transpose()
    |> Nx.reshape({:auto, 4, 2})

  co2_polygon_tensor =
    [
      polygon_tensor,
      xco2 |> Nx.new_axis(-1) |> Nx.tile([4]) |> Nx.new_axis(-1)
    ]
    |> Nx.concatenate(axis: 2)

  case japan_indeces do
    [] ->
      nil
    japan_indeces ->
      Nx.take(co2_polygon_tensor, Nx.tensor(japan_indeces), axis: 0)
  end
end

monthly_polygon_tensor =
  [oco2_dir, target_month_dir]
  |> Path.join()
  |> File.ls!()
  |> Enum.sort()
  |> Enum.map(fn filename ->
    load_japan_polygon_data.(filename)
  end)
  |> Enum.filter(& &1 != nil)
  |> Nx.concatenate(axis: 0)

geometries =
  monthly_polygon_tensor
  |> Nx.to_batched(1)
  |> Enum.map(fn co2_polygon ->
    coordinates =
      co2_polygon[[0, 0..3, 0..1]]
      |> Nx.to_flat_list()
      |> Enum.chunk_every(2)
      |> Enum.map(&List.to_tuple(&1))

    co2 =
      co2_polygon[[0, 0, 2]]
      |> Nx.subtract(min_co2)
      |> Nx.multiply(255 / (max_co2 - min_co2))
      |> Nx.as_type(:u8)
      |> Nx.to_number()

    %Geo.Polygon{
      coordinates: [coordinates],
      properties: %{"xco2" => co2}
    }
  end)
  |> Enum.to_list()

geo_collection = %Geo.GeometryCollection{geometries: geometries}

MapLibre.new(center: {136, 36}, zoom: 3, style: :terrain)
|> MapLibre.add_geo_source("co2_geo", geo_collection)
|> MapLibre.add_layer(
  id: "co2_polygon",
  source: "co2_geo",
  type: :fill,
  paint: [fill_color: jet_expression, fill_opacity: 0.5]
)

富士山の辺りは低く、茨城や千葉の辺りは高そうです

年毎に地図表示する

今度は時系列で変化を見てみましょう

2016年から2022年の毎年12月1日のデータを読み込みます

load_tensors = fn file_path ->
  {:ok, nc_file} = NetCDF.File.open(file_path)

  for var_name <- ["longitude", "latitude", "xco2"], into: %{} do
    tensor =
      nc_file
      |> NetCDF.Variable.load(var_name)
      |> elem(1)
      |> then(&Nx.tensor(&1.value, type: as_nx_type.(&1.type)))

    {var_name, tensor}
  end
end

month_list = ["201612", "201712", "201812", "201912", "202012", "202112", "202212"]

time_series_tensors =
  month_list
  |> Enum.map(fn month_dir ->
    [oco2_dir, month_dir]
    |> Path.join()
    |> File.ls!()
    |> Enum.sort()
    |> Enum.at(0)
    |> then(&[oco2_dir, month_dir, &1])
    |> Path.join()
    |> load_tensors.()
  end)

2016年から2022年の全体で二酸化炭素濃度の最小値、最大値を算出します

full_xco2 =
  time_series_tensors
  |> Enum.map(& &1["xco2"])
  |> Nx.concatenate()

{max_co2, min_co2} =
  {
    full_xco2 |> Nx.reduce_max() |> Nx.to_number(),
    full_xco2 |> Nx.reduce_min() |> Nx.to_number(),
  }

タブ表示で地図にプロットします

time_series_tensors
|> Enum.zip(month_list)
|> Enum.map(fn {tensors, month} ->
  corrected_co2 =
    tensors["xco2"]
    |> Nx.subtract(min_co2)
    |> Nx.multiply(255 / (max_co2 - min_co2))
    |> Nx.as_type(:u8)

  map_source =
    %{
      "latitude" => Nx.to_flat_list(tensors["latitude"]),
      "longitude" => Nx.to_flat_list(tensors["longitude"]),
      "xco2" => Nx.to_flat_list(corrected_co2)
    }

  map =
    MapLibre.new()
    |> MapLibre.add_table_source(
      "xco2",
      map_source,
      {:lng_lat, ["longitude", "latitude"]},
      properties: ["xco2"]
    )
    |> MapLibre.add_layer(
      id: "xco2_layer",
      source: "xco2",
      type: :circle,
      paint: [circle_color: jet_expression]
    )

  {month <> "01", map}
end)
|> Kino.Layout.tabs()

明らかに、毎年二酸化炭素濃度が全体的に上昇しています

二酸化炭素濃度の時系列変化をグラフ化する

グラフでも二酸化炭素濃度の上昇を確認してみましょう

各年の平均二酸化炭素濃度を算出します

mean_xco2_list =
  time_series_tensors
  |> Enum.map(fn tensors ->
    Nx.mean(tensors["xco2"]) |> Nx.to_number()
  end)

グラフ表示します

plot_data =
  mean_xco2_list
  |> Enum.zip(month_list)
  |> Enum.map(fn {mean_xco2, month} ->
    %{
      mean_xco2: mean_xco2,
      year: month |> String.slice(0..3) |> Integer.parse() |> elem(0)
    }
  end)

VegaLite.new(width: 700)
|> VegaLite.data_from_values(plot_data)
|> VegaLite.mark(:line)
|> VegaLite.encode_field(:x, "year", type: :nominal)
|> VegaLite.encode_field(:y, "mean_xco2", type: :quantitative, scale: %{"domain" => [400, 420]})

教科書に出てきそうなぐらい綺麗に比例しています

毎年二酸化炭素濃度は増加しているようです

このことは気象庁のデータからも明らかです

回帰直線を求めて将来の値を予測する

ではこのまま増加すると将来どれくらいまで上がるのか見てみましょう

回帰直線は Evision.fitLine を使って求められます

points =
  mean_xco2_list
  |> Enum.zip(month_list)
  |> Enum.map(fn {mean_xco2, month} ->
    [month |> String.slice(0..3) |> Integer.parse() |> elem(0), mean_xco2]
  end)
  |> Nx.tensor()

[vx, vy, x, y] =
  points
  |> Evision.fitLine(Evision.Constant.cv_DIST_L2(), 0, 0.01, 0.01)
  |> Evision.Mat.to_nx(EXLA.Backend)
  |> Nx.to_flat_list()

slope = vy / vx
intercept = y - (slope * x)

{slope, intercept}

実行結果は以下のようになります

{2.3805762028504414, -4396.929740273791}

• 傾き: 2.3805762028504414
• 切片: -4396.929740273791

つまり、西暦と二酸化炭素濃度は以下のような式で表されます

$ XCO_2 = 2.3805762028504414Year - 4396.929740273791 $

このまま行くと西暦2100年にはどうなるかというと

xco2_predictions =
  2023..2100
  |> Enum.map(fn year ->
    %{
      mean_xco2: slope * year + intercept,
      year: year
    }
  end)

xco2_predictions
|> Enum.at(-1)

実行結果

%{mean_xco2: 602.280285712136, year: 2100}

600ppm を超えるようです

ではこれをグラフに表してみましょう

predictions_data = plot_data ++ xco2_predictions

VegaLite.new(width: 700)
|> VegaLite.data_from_values(predictions_data)
|> VegaLite.mark(:line)
|> VegaLite.encode_field(:x, "year", type: :nominal)
|> VegaLite.encode_field(:y, "mean_xco2", type: :quantitative)

結構な増加率であることが分かると思います

二酸化炭素濃度が上昇するとどの程度温暖化するのか、まではこのデータだけでは分かりませんが

まとめ

改めて自分で二酸化炭素濃度を可視化することで、はっきりとその変化を実感することができました

こういったデータサイエンスを学校の授業にも取り入れられるといいですね