熊本県の植物の分布を調べてみる

# はじめに
今回は、Pythonを利用して、ANNIR-2によって撮影された熊本県の画像を取得し、植物の分布を調べてみました。初心者なので至らない点も多いと思いますが、温かく見守って頂くと幸いです

# 衛生データを所得
### 解析に必要ライブラリをインポートとAPIトークンを入力
```python
import os
os.environ['USE_PYGEOS'] = '0'
import requests
import json
import glob
import rasterio.mask
import fiona
import cv2

import numpy as np
import rasterio as rio
import matplotlib.pyplot as plt
import geopandas as gpd
from rasterio.plot import show
from rasterio.plot import show_hist
```
### 次にTellusのAPIトークンを入力します
```python
token = '#################'
```

### 2011年4月6日にAVNIR-2によって撮影された熊本県の画像をとってくる
Tellus Travelarで必要な衛生データに関する情報を取得し、それを入力する

```python
dataset_id = 'ea71ef6e-9569-49fc-be16-ba98d876fb73'

scenes = [
    '6a4d2a0b-31a4-428b-af0b-b64c7489e219'
]
```
### データセットIDとシーンIDから、そのシーンに含まれるファイル名を取得する関数をかく

```python
# データセットIDとシーンIDから、そのシーンに含まれるファイル名を取得する
def get_dataset_data_by_id_files(base_path, token, dataset_id, data_id):
    print(base_path, token, dataset_id, data_id)
    url = '{}/api/traveler/v1/datasets/{}/data/{}/files/'.format(
        base_path, dataset_id, data_id)
    headers = {
        'Authorization': 'Bearer ' + token,
        'content-type': 'application/json'
    }
    r = requests.get(url, headers=headers)
    assert r.status_code == 200
    return json.loads(r.content)
```

### ファイルIDを指定して、ダウンロードURLを取得する関数をかく
```python
# ファイルIDを指定して、ダウンロードURLを取得する
def get_dataset_data_by_id_files_by_id_download_url(base_path, token, dataset_id, data_id, file_id):
    url = '{}/api/traveler/v1/datasets/{}/data/{}/files/{}/download-url/'.format(
        base_path, dataset_id, data_id, file_id)
    headers = {
        'Authorization': 'Bearer ' + token,
        'content-type': 'application/json'
    }
    r = requests.post(url, headers=headers)
    assert r.status_code == 200
    return json.loads(r.content)
```
### 衛生データのダウンロードを行う関数をかく

```python
# 衛生データのダウンロードを行う
def download(scenes, dataset_id, token, dist='./', base_path='https://www.tellusxdp.com'):
    for scene_id in scenes:
        # データセットIDとシーンIDから、そのシーンに含まれるファイル名を取得する
        files = get_dataset_data_by_id_files(
            base_path, token, dataset_id, scene_id
        )
        print(files)
        
        rawdata = files['results']
        path = os.path.join(dist, scene_id)
        
    if len(rawdata) > 0:
        for file in rawdata:
            file_id = file['id']
            file_name = file['name']
            
            file_path = os.path.join(path, file_name)
            
            # ファイルIDを指定して、ダウンロードURLを取得する
            download_url = get_dataset_data_by_id_files_by_id_download_url(
                base_path, token, dataset_id, scene_id, file_id)['download_url']
            
            r = requests.get(download_url, stream = True)
            
            if not os.path.exists(path):
                os.makedirs(path)
                
            with open(file_path, 'wb') as f:
                for chunk in r.iter_content(chunk_size=1024):
                    if chunk:
                        
                        f.write(chunk)
                        f.flush()
```

### データをダウンロード
```python
# データDL
download(scenes, dataset_id, token)
```

# NDVIを算出し、NDVI画像を出力

### NDVIを算出する関数をかく
```python
# NDVIの算出
def mk_ndvi_rio(red_img, nir_img, red_gain, red_offset, nir_gain, nir_offset, out_ndvi):
    red = rio.open(red_img)
    nir = rio.open(nir_img)
    
    red_arr = red.read().astype(float)
    nir_arr = nir.read().astype(float)
    
    kwargs = {
        "width": red.width,
        "height": red.height,
        "count": 1,
        "crs": red.crs,
        "transform": red.transform,
        "dtype": rio.float32
    }
    
    ## DN2輝度変換
    red_arr = red_arr * red_gain + red_offset
    nir_arr = nir_arr * nir_gain + nir_offset
    
    # NDVI算出
    ndvi = np.where(((red_arr != 0) & (nir_arr != 0)), ((nir_arr - red_arr) / (nir_arr + red_arr)), 0)
    
    write_ndvi = rio.open(out_ndvi, "w", driver="Gtiff", **kwargs)
    write_ndvi.write(ndvi)
    write_ndvi.close()
    
    return(out_ndvi)
```

### NDVI画像を作成
```python
dir_path = "D:\\そら畑\\植物の分布から田畑の利用状況を推測する衛星データの利活用方法\\6a4d2a0b-31a4-428b-af0b-b64c7489e219"


subdir_list = os.listdir(dir_path)
red_path = glob.glob(dir_path + "/IMG-03-*.tif")
nir_path = glob.glob(dir_path + "/IMG-04-*.tif")
hdr_path = glob.glob(dir_path + "/HDR-*.txt")
hdr_path = str(hdr_path[0])

# HDR情報から、ゲイン、オフセット情報を抽出
hdr_f = open(hdr_path,"r")
hdr_data = hdr_f.read()
hdr_f.close

# gain, offset
red_gain = float(hdr_data[1753:1761])
red_offset = float(hdr_data[1761:1769])
nir_gain = float(hdr_data[1769:1777])
nir_offset = float(hdr_data[1777:1785])

red_path = red_path[0]
nir_path = nir_path[0]
print(red_path, nir_path)

#NDVI画像の作成
base_name = os.path.basename(red_path)
ndvi_file = "NDVI-" + base_name[7:]
ndvi_path = os.path.join(dir_path, ndvi_file)
ndvi_path

mk_ndvi_rio(red_path, nir_path, red_gain, red_offset, nir_gain, nir_offset, ndvi_path)
```

### NDVI画像とカラー画像を出力
```python
ndvi_20110406 = rio.open("D:\\そら畑\\植物の分布から田畑の利用状況を推測する衛星データの利活用方法\\6a4d2a0b-31a4-428b-af0b-b64c7489e219\\NDVI-ALAV2A276862940-OORIRFU-D083P3-20110406-002.tif")
RGB_20110406 = cv2.imread('ALAV2A276862940-OORIRFU-D083P3-20110406-002_thumb.png')
ndvi_20110406 = ndvi_20110406.read(1)

plt.figure(figsize=(10, 10))

plt.subplot(1, 2, 1)
plt.imshow(cv2.cvtColor(RGB_20110406, cv2.COLOR_BGR2RGB))

plt.subplot(1, 2, 2)
plt.imshow(ndvi_20110406, cmap='RdYlGn', vmin=-1, vmax=1)

plt.show()

type(RGB_20110406)
```

出力
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3494144/34f634f7-57e4-d805-634b-bf0db7389c8c.png)

左がカラー画像で右がNDVI画像です。カラー画像で植物の覆われているところは、NDVI画像では緑で表示されています。逆に植物が活発じゃないところは、オレンジで表示されています。

# まとめ
今回は、熊本県の植物の分布を調べてみました。これからもいろんな衛生データ解析をやっていきたいと思います。

# 参考文献及び引用
植物の分布から田畑の利用状況を推測する衛星データの利活用方法: https://sorabatake.jp/28126/