Google Earth Engineで道路周辺の都市化（建物）の変化を見る

# 1. はじめに
衛星データを活用して、ネパールのシンズリ道路周辺の建物の変化を見たいと思います。

シンズリ道路については、[こちら](https://qiita.com/ishimotojuri/items/0d14be0ada5a9cf2a955)をご覧ください。

# 2. データ
土地利用データから建物を抽出します。

土地利用データには、[Dynamic World V1](https://developers.google.com/earth-engine/datasets/catalog/GOOGLE_DYNAMICWORLD_V1)を使用します。
Dynamic World は、9つ​​の土地利用の書類と、その確率、ラベル情報を含む、解像度10mのデータセットです。

# 3. 分析手順
## 3.1. 関心領域の設定
関心領域（シンズリ道路）のシェープファイルを用意します。
シェープファイルの作成方法は[こちら](https://qiita.com/ishimotojuri/items/a74ecbde027ca3c06fe7)をご覧ください。

関心領域のシェープファイルをGoogle Eart Engine（GEE）にインポートします。
インポートの方法は[こちら](https://developers.google.com/earth-engine/guides/table_upload)をご参照ください。

<INS>シンズリ道路周辺のポリゴンを作成します</INS>。
マップ上の左にある、「長方形を描画」をクリックして、適当にポリゴンを作成します。
シンズリ道路周辺の分析には、このポリゴンを使用します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/e84045dd-58f9-d104-dde0-4348603fa22d.png)

## 3.2. 土地利用データのインポート
<INS>GEE上で2016年と2023年の土地利用データを取得します</INS>。
```js
// Import the Satellite data
// 2016
var pre_data = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1')
             .filterDate('2016-01-01', '2016-12-31')
             .filterBounds(geometry)
             .select('label')
             .reduce(ee.Reducer.mode()); 

// 2023
var post_data = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1')
             .filterDate('2023-01-01', '2023-12-31')
             .filterBounds(geometry)
             .select('label')
             .reduce(ee.Reducer.mode()); 
```

## 3.3. 変化の可視化
<INS>インポートした画像から、建物と分類されたピクセルを抽出します</INS>。      
```js
// Extract builtAreas
var pre_builtArea = pre_data.eq(6);
var post_builtArea = post_data.eq(6);
```
<INS>2023年と2016年の差分を計算し、変化をマップに表示させます</INS>。   
```js
// Calculate the change
var diff_builtArea = post_builtArea.subtract(pre_builtArea);

var loss_builtArea = diff_builtArea.eq(-1); // extract loss area (2016-2023)
var gain_builtArea = diff_builtArea.eq(1); // extract gain area (2016-2023)

// Visualize the change
var loss_builtAreaMasked = loss_builtArea.selfMask();
var gain_builtAreaMasked = gain_builtArea.selfMask();

Map.addLayer(loss_builtAreaMasked.clip(geometryBuffer),
{palette: ['blue']}, 
'[blue]: Loss of buitAreas');

Map.addLayer(gain_builtAreaMasked.clip(geometryBuffer),
{palette: ['red']}, 
'[red]: Gain of buitAreas');

// Display relevant geometries on the map
Map.addLayer(bphighway, {'color': 'black'}, '[black]: bphighway');
Map.centerObject(bphighway,10);
```        

青は建物が減少（2016年に建物があったが、2023年にない）、赤は建物が増加（2016年に建物がなかったが、2023にある）していることを示しています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/91165349-edf9-849f-e800-71678a3d8b4b.png)

## 3.4. 変化の画像のエクスポート
<INS>変化の画像をエクスポートします</INS>。
```js
// Export the image
var diff_builtArea_clip = diff_builtArea.clip(geometry);

Export.image.toDrive({
  image: diff_builtArea,
  description: 'diff_builtArea',
  folder: 'Nepal', 
  scale: 10,
  region: geometry,
  maxPixels: 1e12
});
```

上記のコードを実行すると、右側の「Tasks」画面で各ファイル名が表示されます。
ファイル名の「Run」ボタンをクリックすると、画面が表示されます。

その画面の設定条件を確認し、下の「Run」ボタンをクリックするとエクスポートが開始されます。

ダウンロードしたデータはGoogleドライブの任意のフォルダに保存されます。
 
# 4. 分析結果
ダウンロードした画像は、下図のように、QGIS上で加工できます。
シンズリ道路沿いの町のシェープファイルも追加しています。

図から、シンズリ道路周辺で建物が増えていることが確認されました。さらに、シンズリ道路周辺のうち特にドリケル、バクンベデシ、シンズリバザール、バルディバスという町の周辺で建物が増えていることがわかりました。
![5_builtArea_change_road.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/e092fde5-4b1a-edfd-5990-493548f2003b.png)


# 5. おわりに
シンズリ道路周辺において、2016年と比較して2023年に建物が増加していることが明らかになりました。

現地調査では、シンズリ道路開通後、周辺地域から沿道沿いの町に人々が来るようになり、新しい商店やビジネスが増加したといったことを聞きました。このような定性的な情報が衛星データによって補完されたと考えられます。