森林面積の算出

# はじめに
この記事では[MODIS Land Cover Type Yearly Global 500m](https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MCD12Q1)を使用して、森林面積を算出する方法を紹介します。
このデータは土地被覆図です。データ期間は2001年1月1日～2020年1月1日までで、解像度は500mです。ここでは、International Geosphere-Biosphere Programme (IGBP)に基づく森林5クラスの総面積を集計することで、森林面積を算出します。

# 分析環境
Google Earth Engineを使用しました。


# 手順
1\. 関心領域のインポート（今回はマラウイとします）
2\. MODISデータのインポート
3\. MODISデータの可視化
4\. 2001年と2019年の森林面積の算出
5\. 2001年と2019年の変化を可視化


# 1\. 関心領域のインポート
```js:sample.js
var mi = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
         .filter(ee.Filter.eq("country_co", "MI"));

Map.centerObject(mi, 6);
Map.addLayer(mi, {}, "Malawi");
```

<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/fab85eea-4f8a-f259-933d-2ce405b2bab8.png">

# 2\. MODISデータのインポート
```js:sample.js
var modis_01 = ee.ImageCollection("MODIS/006/MCD12Q1")
            .filterDate("2001-01-01", "2001-12-31")
            .select("LC_Type1")
            .median()
            .clip(mi);



var modis_19 = ee.ImageCollection("MODIS/006/MCD12Q1")
            .filterDate("2019-01-01", "2019-12-31")
            .select("LC_Type1")
            .median()
            .clip(mi);
```

# 3\. MODISデータの可視化
```js:sample.js
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_01, igbpLandCoverVis, 'IGBP Land Cover Malawi in 2001');
Map.addLayer(modis_19, igbpLandCoverVis, 'IGBP Land Cover Malawi in 2019');
```

2001年と2019年のデータを可視化してみました。図を比べると、中央より上の森林が少し減少しているように見えます。

<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/64d52862-9ef7-168d-202f-d2a3d960f06d.png">
<img width="296" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/6a3052f9-ac2a-bfb3-9989-5bdf0828c17f.png">

左：2001年のデータ、右：2019年のデータ

# 4\. 2001年と2019年の森林面積の算出
```js:sample.js
var tree_cover_mask_01 = modis_01.lt(6);
var tree_cover_01 = modis_01.updateMask(tree_cover_mask_01);

Map.addLayer(tree_cover_01, igbpLandCoverVis, "Filtered Forest Map Malawi in 2001");

var tree_cover_mask_19 = modis_19.lt(6);
var tree_cover_19 = modis_19.updateMask(tree_cover_mask_19);

Map.addLayer(tree_cover_19, igbpLandCoverVis, "Filtered Forest Map Malawi in 2019");
```

<img width="250" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/0db4f207-bca0-0485-fb44-3ccc93d8d4f4.png">
<img width="239" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/743924be-7b8b-6ea5-5033-bbda1a3a1290.png">

左：2001年のデータ、右：2019年のデータ

森林のピクセル数を算出します。
```js:sample.js
var forest_pixels_01 = tree_cover_01.reduceRegion({
  reducer:ee.Reducer.sum(),
  geometry: mi.geometry(),
  scale:500,
  maxPixels: 1e10
});

print('Total of pixcel of forest in Malawi 2001',forest_pixels_01);

var forest_pixels_19 = tree_cover_19.reduceRegion({
  reducer:ee.Reducer.sum(),
  geometry: mi.geometry(),
  scale:500,
  maxPixels: 1e10
});

print('Total of pixcel of forest in Malawi 2019',forest_pixels_19);
```

右のコンソールにそれぞれ以下の数字が表示されました。

2001年のピクセル数：111607.662745098
2019年のピクセル数：100930.88627450979

解像度は500メートル四方なので、上記の数字に250000mをかけ、面積を算出します。
そうすると、

2001年の森林面積(k㎡)：27,901
2019年の森林面積(k㎡)：25,232

となりました。よって、森林が減少していることが計算からも明らかになりました。

# 5. 2001年と2019年の変化を可視化
2001年に森林があった箇所で、2019年になければ（森林損失）、赤に表示させます。

```js:sample.js
var change = tree_cover_01.subtract(tree_cover_19);
var change_mask = change.gt(0.1);
var visParams_loss = {palette:["red"]};

Map.addLayer(change.updateMask(change_mask), visParams_loss, "Forest loss in Malawi from 2001 to 2019");
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/f0185455-3003-23b9-c44a-eaf7cd5ed552.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/eccdbdf9-f693-500a-9754-1e15db662e8f.png)

上：マラウイ全体図
下：中央部分の拡大図

上記で示したコードは[こちら](https://code.earthengine.google.com/d832e3454b018dc290adc7bf3aa9fd6c)から見れます。

# おわりに
上記の通り、森林面積の算出及び森林が損失した箇所の可視化により、2001年から2019年にかけて一部の森林が損失したことがわかりました。ただし、今回使用したのは解像度が粗いデータ（500m）ですので、現実との乖離があるかもしれない点には注意が必要です。今後、より解像度の高いデータを使って、森林面積の変化を捉えてみたいと思います。