はじめに

コンゴ民主共和国（以下、コンゴ民）は、国別ではブラジルに次いで世界で二番目に多い熱帯雨林を有していますが、違法伐採を含む商業伐採や農地開発、鉱山開発等により毎年多くの森林が失われています。

米国の世界資源研究所がGoogleなどの企業や大学、研究所、国際機関などと連携して作成した「グローバル・フォレスト・ウォッチ」は、人工衛星などで分析した最新の森林減少面積を無料で公開しています。これによると、2020年から2021年のコンゴの熱帯原生林の減少面積は499khaでした。これは東京都（約22万ha）の約2.3倍。たったの一年で、「地球の肺」と呼ばれるかけがえのない森が、大量に消えてしまっています。

こうした中、コンゴ民の森林を守ろうと、中部アフリカ森林イニシアティブ（Central Africa Forest Initiative: CAFI）が立ち上げられ、2016年～2020年の期間にコンゴ民で約2億ドル（約270億円）の資金コミットメントがなされました。さらに今後、約5億USD（約680億円）の資金コミットメントがなされる予定です。このイニシアティブの枠組みの中で、多くの援助機関が支援を行っています。

日本のODA実施機関であるJICAも、2019年4月～2024年4月までの5年間の計画で、「国家森林モニタリングシステム運用・REDD＋パイロットプロジェクト」を実施し、成果１：国家森林モニタリングシステム（NFMS）の運用能力強化、成果２：パイロット州（クウィル州）におけるREDD＋事業の実施（「PIREDD Kwilu」）の2つの成果を設定して取り組んでいます。このプロジェクトでは、CAFI資金の一部（約4億円）を受託し、同資金も活用しながら、クウィル州の約250か村を対象にアグロフォレストリーや森林保全活動に取り組んでいます。

このように多くの機関が支援を行っており、プロジェクトの森林減少抑制への効果を検証することは非常に重要だと思われます。
今回は上記のJICAのプロジェクトの中間時点での効果を検証してみたいと思います。

分析環境

Google Earth Engineを使用しました。

手順

1.　関心領域のインポート（今回はパイロットプロジェクトを行っているクウィル州を対象とします）
2.　2時点の森林面積の算出
3.　プロジェクト対象地域の森林面積の変化の算出
4.　プロジェクト非対象地域の森林面積の変化の算出
5.　プロジェクト対象地域と非対象地域の森林面積の変化の差の算出

なお、効果検証の場合はボリュームが多くなるため、今回は複数に分けて記事を執筆・投稿したいと思います。
一番目の投稿は1～2までです。

1. 関心領域のインポート

1.1.　行政区境界データの取得

クウィル州の領域を示すデータ（行政区境界データ）を収集します。Google Earth Engineで使えるデータは2015年時点のものなのですが、コンゴ民の行政区はそれ以降変わっているようです。2015年のデータには、お目当てのクウィル州が見当たりませんでした。そのため、以下の別のウェブサイトから入手します。

The Humanitarian Data Exchange（HDX）

Humanitarian Data Exchange（HDX）は、国際連合人道問題調整事務所（OCHA）が運営するオープンデータプラットフォームです。HDXは、人道支援や災害対応に関する情報やデータを共有し、人道支援の効果を高めることを目的としています。

このサイトでは、コンゴ民のDistrictレベル（level 2）までのデータが取得できます。
プロジェクトはクウィル州を対象としていますが、2019年に開始して以降、厳密にはクウィル州を構成する5つのDistrict（Bagata, Masi-manimba, Bulungu, Idiofa, Gungu）のうちGunguを除く4つのDistrictで介入を行っています（Gunguへの介入は今後予定）。あとで介入地域と非介入地域での変化を比較することができるよう、Districtレベルのデータ（シェープファイル形式）をダウンロードします。

QGISを立ち上げ、QGISの画面にシェープファイルをドラッグ・ドロップして、レイヤを右クリック⇒属性テーブルからクウィル州が入っているかを確認。下図の通り、クウィル州が入っていることが確認できました。

1.2.　GEEへの行政区境界データのインポート

GEEのコードエディターを開き、左ウィンドウのAssets⇒NEW⇒Table Upload⇒Shapefilesをクリックし、先ほど確認したAdm2（District）のファイル（拡張子がshp, dbf, shx）をドラック・ドロップ。Asset IDはわかりやすい名前にしておきましょう。その他は基本はデフォルトの設定のままで問題ありません。
UPLOADを押して完了です。何分か待つことになります。実行状況は右ウィンドウのTasksで確認できます。

UPLOADが完了したら、左ウィンドウのAssets⇒LEGASY ASSESTSで先ほど設定したIDがあるはずです（でていなかったらNEWの右隣のRefreshボタンを押してください）。
IDからシェープファイルのTable IDや属性が確認できます。このTable IDを呼び込むことで、Google Earth Engineで使うことができます。

コードは以下の通りです。FeatureCollection(Table ID).filter(ee.Filter.eq("使いたい行政区域がある列の名前", "使いたい行政区域"))と設定してください。

sample.js

var Kwilu = ee.FeatureCollection("users/ishimotojuri/DRC_adm2")
         .filter(ee.Filter.eq("ADM1_FR", "Kwilu"));

Map.centerObject(Kwilu, 6);
Map.addLayer(Kwilu, {}, "Kwilu");

2. 2時点の森林面積の算出

森林面積の算出には、以下の方法があります。

森林境界や土地利用区分に関する地理情報システム（GIS）データを組み合わせて、森林境界の面積を算出する。
衛星画像から作成された土地被覆マップを活用し、森林面積を算出する。（他人が作ったものを活用する）
衛星画像で機械学習を行い、土地被覆マップを作成し、森林面積を算出する。（自分で作成する）

今回は2の方法で森林面積を算出したいと思います。
衛星データは各年の森林面積を示すMODIS Land Cover Type Yearly Global 500mを使用します。
このデータは、MODIS(Moderate Resolution Imaging Spectroradiometer)センサーを使用して撮影された地球の全地表面について、年次で500mの空間解像度で提供される地表被覆タイプのデータセットです。このデータセットは、世界の陸域を21の地表被覆タイプに分類しています。これらのタイプには、森林、草原、サバンナ、湿地、耕地、都市、砂漠などが含まれます。
Google Earth Engineでは、2001年～2020年までの情報が取得可能で、解像度は500mです。森林はさらに5クラスに分かれており、これらの総面積を集計することで「森林面積」を算出することが可能になっています。

2.1.　MODIS土地被覆マップでの森林面積の算出

2.1.1.　プロジェクト開始前（2018年）の森林面積の算出

プロジェクトが開始する前の2018年の森林面積を算出します。まずはクウィル州全体で土地被覆図を表示してみます。

sample.js

var modis_18 = ee.ImageCollection("MODIS/061/MCD12Q1")
                 .filterDate("2018-01-01", "2018-12-31")
                 .select("LC_Type1")
                 .median()
                 .clip(Kwilu);

var igbpLandCoverVis = {
  min: 1.0,
  max: 17.0,
  palette: [
    '05450a', '086a10', '54a708', '78d203', '009900', 'c6b044', 'dcd159',
    'dade48', 'fbff13', 'b6ff05', '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c',
    '69fff8', 'f9ffa4', '1c0dff'
  ],
};

Map.addLayer(modis_18, igbpLandCoverVis, 'MODIS Land Cover in Kwilu in 2018');

黄色はサバンナ（低木が植えられている）、黄緑色は草原、緑色は森林なので、全体的にサバンナ地域が広がっていることがわかります。実際に現場に行った際も、昔は森があったが、伐採してしまい今はサバンナ地域（乾燥して低木がまばらにある）というところをいくつか見かけました。

次に森林と分類されているセルだけ抽出して表示させます。これには、updateMaskという関数を使います。
updateMaskは、与えられたマスク画像が1（有効なピクセル）の場合、元の画像のピクセル値を保持します。マスク画像が0（無効なピクセル）の場合、元の画像のピクセル値はマスク処理され、0に設定されます。例えば、森林面積を計算するために、森林以外（ここではクラスの6未満）の部分を0（マスク処理）にし、森林部分（ここではクラスが1~5）のピクセル値のみを使用する場合、updateMaskを使用して、森林以外の部分をマスク処理（特定の部分を切り取る、排除する）することができます。

例えば、以下のコマンドでは、1行目でマスク画像を作成しており、クラスの6未満（1～5の森林地域）をもつピクセルの値を1、それ以外のクラスを持つピクセルの値を0に設定します。
次に2行目でupdateMaskの関数を使用し、マスク画像が1の場合に元の画像（modis_18）のピクセル値はそのまま（森林地域はそのまま）、マスク画像が0の場合に元の画像のピクセルは排除（森林地域以外のピクセル値は0になる）させます。

sample.js

var tree_cover_mask_18 = modis_18.lt(6);
var tree_cover_18 = modis_18.updateMask(tree_cover_mask_18);

Map.addLayer(tree_cover_18, igbpLandCoverVis, "Filtered Forest Map in Kwilu in 2018");

すると、画像は以下のようになります。森林地域だけが緑として残りました。（境界のレイヤーも透過性を変えて残してあります）

森林のピクセルを抽出したら、森林面積を算出します。面積の算出にはee.Image.pixelArea()というメソッドを使います。
これを使用すると、各ピクセルの面積が正確に計算され、画像の処理（森林面積の計算や土地利用変化の解析）を行うことができるようになります。
単位は平方メートル（m²）です。また、面積の計算には、画像の投影法や解像度などの情報が使用されます。

sample.js

var forest_area_m2_18 = tree_cover_18.multiply(ee.Image.pixelArea());

var total_forest_area_m2_18 = forest_area_m2_18.reduceRegion({
  reducer:ee.Reducer.sum(),
  geometry:Kwilu.geometry(),
  scale:500,
  maxPixels: 1e10
});

print('Total of forest area m2 in Kwilu in 2018;', total_forest_area_m2_18);
//11391628788.319115m2
//11,391km2
//1,139,100ha

2.1.2.　プロジェクト開始後（2021年）の森林面積の算出

上記と同様の手順でプロジェクト開始後（2021年）の森林面積を算出します。
コードは以下の通りとなります。

sample.js

var modis_21 = ee.ImageCollection("MODIS/061/MCD12Q1")
                 .filterDate("2021-01-01", "2021-12-31")
                 .select("LC_Type1")
                 .median()
                 .clip(Kwilu);
                 

Map.addLayer(modis_21, igbpLandCoverVis, 'MODIS Land Cover in Kwilu in 2021');


var tree_cover_mask_21 = modis_21.lt(6);
var tree_cover_21 = modis_21.updateMask(tree_cover_mask_21);

Map.addLayer(tree_cover_21, igbpLandCoverVis, "Filtered Forest Map in Kwilu in 2021");


var forest_area_m2_21 = tree_cover_21.multiply(ee.Image.pixelArea());

var total_forest_area_m2_21 = forest_area_m2_21.reduceRegion({
  reducer:ee.Reducer.sum(),
  geometry:Kwilu.geometry(),
  scale:500,
  maxPixels: 1e10
});

print('Total of forest area m2 in Kwilu in 2021;', total_forest_area_m2_21);
//11202735577.149141
//11,202km2
//1,120,200ha

2.1.3.　注意事項

なお、面積を算出する際に気をつけなければならないのが、画像をメートルの単位にすることです。メートルなどの距離を算出するには、UTM（Universal Transverse Mercator coordinate system）座標にする必要があります。UTM座標では、座標は東西方向（x座標）と南北方向（y座標）で表されます。x座標は、地域の中心のグリッド線からの距離で表され、y座標は赤道からの距離で表されます。UTM座標系は、単位がメートルです。

具体的には、以下の手順で画像をメートル単位にすることができます。

1. ee.Image.projection()メソッドを使用して、画像の投影情報を取得する。
2. ee.Image.reproject()メソッドを使用して、画像の投影法を変更する。

結果は以下の通りとなりました。
変更前：11,391km2、変更後（UTM座標）：11,368km2

投影法はキンシャサ周辺のUTM、EPSG:32733にしました。座標を変えて計算し直しましたが、上記とさほど変わらない結果になりました。
これは、Google Earth Engineが元の座標系であるWGS84（EPSG:4326）とUTM座標系の間で自動的に変換を行うためだそうです。
この変換は非常に正確であるため、変換前後の座標の位置にほとんど差異がない場合があります。
ただし、広範囲にわたる画像や地理情報を処理する場合には、座標系の変換による精度の損失が考慮される必要があります。特に、測地系が異なる場合や、座標系の歪みが大きい場合には、座標系の変換が精度に影響を与える可能性があります。そのため、広範囲にわたる画像や地理情報を処理する際には、座標系の選択や変換方法について慎重に考慮する必要があります。

sample.js

// Get the projection information
var projection = tree_cover_18.projection();

// Get the CRS of the image
var crs = projection.crs();

// Print the CRS of the image
print('CRS:', crs);

var utmProjection = ee.Projection('EPSG:32733');

// Reproject the image to UTM projection
var tree_cover_18_UTM = tree_cover_18.reproject({
  crs: utmProjection,
  scale: 500,
});

print('UTM座標系に再投影された画像', tree_cover_18_UTM);

var forest_area_m2_18_UTM = tree_cover_18_UTM.multiply(ee.Image.pixelArea());

var total_forest_area_m2_18_UTM = forest_area_m2_18_UTM.reduceRegion({
  reducer:ee.Reducer.sum(),
  geometry:Kwilu.geometry(),
  scale:500,
  maxPixels: 1e10
});

print('Total of forest area m2 in Kwilu in 2018;', total_forest_area_m2_18_UTM);
//11368909219.30993
//11,368km2
//1,136,800ha

おわりに

クウィル州全体の森林面積を算出してみると、2018年（プロジェクト開始前）は11,391km2、2021年（プロジェクト開始後）は11,202km2であり、全体的に減少傾向にあることがわかりました。
ただし、プロジェクトは州全体に介入しているわけではないため、プロジェクトの効果を厳密に検証できていない可能性があります。
そのため、次回はプロジェクトの介入地域と非介入地域を特定し、その差を検出するインパクト評価を行ってみたいと思います。

また、今回使用したModisのデータセットでは、2021年までのものしか入手できなかったため、Sentinel-2で機械学習を行い、2021年以降の森林面積を算出したいと考えています。