Pythonで衛星画像を触ってみる その２

続き

前回までに衛星画像をローカル環境にダウンロードしてきました。本記事ではダウンロードした衛星画像に対して、GDALを利用して何かしらの操作をしていきたいと思います。

参考サイト

環境・前提条件

• その１の記事の内容に沿って、ローカル環境にDLした画像がある。
• JupyterNotebookを利用する
• Gdalがインストールされていること
  　→　参考：https://info.qchizu.xyz/knowledge/knowledge-gdal-install/

やってみる

必要なライブラリをインポートします。

# import necessary libraries
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

ダウンロードした画像の一部を切り出す関数を作成します。

#　Geotiffから指定範囲（緯度経度）を切り出す関数
def cut4_geotiff(IMG_IN, IMG_OUT, minX, maxX, minY, maxY):
    gdal.UseExceptions()
    gdal.Translate(IMG_OUT, IMG_IN, projWin=[minX,maxY,maxX,minY])
    return print("Extraction process was completed! : ", IMG_IN)

ダウンロードしたファイルをそれぞれ指定して、切り出した上で、赤青合成を行い、PNGで表示します。このコードの中で上記で作成した切り出し用の関数を指定しています。

# ダウンロードしたファイルを指定して、切り出した上で、赤青合成を行い、PNGで表示
if __name__ == '__main__':
# 処理対象ファイルの指定
    IMG_PATH_1 = "./AS200355208498-180908_webcog.tif"
    IMG_PATH_2 = "./AS201462808498-200905_webcog.tif"
    
    # AOI抽出後(切り出し後)のGeotiffの保存先の指定
    IMG_PATH_1_CUT = "./out/AS200355208498-180908_webcog_cut.tif"
    IMG_PATH_2_CUT = "./out/AS201462808498-200905_webcog_cut.tif"
    IMG_JPG = "./out/out.jpeg"
    
    # AOI(切り出し部分)の四隅緯度経度指定
    minX = 141.93507323883878
    minY = 42.75282111399301
    maxX = 141.94639048725506
    maxY = 42.75848804932892
    
    # AOIの抽出処理（四隅）: 切り出し処理
    cut4_geotiff(IMG_PATH_1, IMG_PATH_1_CUT, minX, maxX, minY, maxY)
    cut4_geotiff(IMG_PATH_2, IMG_PATH_2_CUT, minX, maxX, minY, maxY)
    
    # IMG_1読込
    ds_1 = gdal.Open(IMG_PATH_1_CUT, gdal.GA_ReadOnly)
    image_1 = np.array([ds_1.GetRasterBand(i + 1).ReadAsArray() for i in range(ds_1.RasterCount)])
    del ds_1 #不要メモリ解放

    # IMG_2読込
    ds_2 = gdal.Open(IMG_PATH_2_CUT, gdal.GA_ReadOnly)
    image_2 = np.array([ds_2.GetRasterBand(i + 1).ReadAsArray() for i in range(ds_2.RasterCount)])
    del ds_2 #不要メモリ解放

    # 画像間のピクセル数を揃える処理(今回は必要なし)
    x_min = min([image_1.shape[1], image_2.shape[1]]) 
    y_min = min([image_1.shape[2], image_2.shape[2]]) 
    # print("MinX : ", image_1.shape[1], image_2.shape[1])
    # print("MinY : ", image_1.shape[2], image_2.shape[2])

    # 出力用配列の準備
    image = np.zeros((x_min,y_min,3)).astype(np.uint8)
    image[:,:,0] = image_2[0, :x_min, :y_min]
    image[:,:,1] = image_1[0, :x_min, :y_min]
    image[:,:,2] = image_1[0, :x_min, :y_min]
    # print(image_1.shape, image_2.shape, image.shape)

    # 不要メモリの解放
    del image_1
    del image_2

    # 出力画像の保存
    pil_img = Image.fromarray(image)
    # print(pil_img.mode)
    pil_img.save(IMG_JPG)

    print("DONE!")

この関数実行時に、私の環境ではエラーが発生し、記事作成時点ではまだ解決ができていません。
エラーについてのQAをTeratailでしてますので、こちらをウォッチしつつ、解決したら本記事も更新しようと思っています。
https://teratail.com/questions/hi5rrkeq01dvft?_from_e=e_qctq_showd

処理が正常に完了されると、以下のような画像が出力されます。

引用：https://sorabatake.jp/25931/

発災直後2018-09-08（青）と発災から2年後2020-09-05（赤）で比較したASNARO-2画像
Credit : NEC Corpotration

以上、ここまでがダウンロードした画像に対する捜査でした。

おわりに

TellusAPIを利用して衛星画像を取得してみて、衛星画像を用いた解析などをPythonで行うには、解析の知識（MLなど）に加えて、衛星画像そのものに対する知識も必要だと感じました。また、衛星画像は（とくに解像度が高い画像は）非常にサイズも大きいので、そのあたりをどのようにハンドルするか（グレースケールにする、必要な部分を切り出しする、など）も考慮する必要があることがわかりました。衛星画像と聞くと「人工衛星」＝「宇宙」＝「ロケットで打ち上げ」など、色々とワクワクするような用語が思い浮かんだりするので、引き続き学習を重ねていこうと思います。