本稿は、pythonを使用して衛星画像をどのように処理するかについてまとめます。

①使用するライブラリ
numpyとgdalとmatplotlib
※numpyは計算処理のため
※gdalは衛星画像用のライブラリ。ほかにもrasterというものがあるが、私の作業場の親和性はgdalのほうが高かった。

gdalのインストールの仕方について、まとめた記事があるのでそちらを参照してください。

②作業
手順としては、
gdalで衛星画像を読み込み
→numpyで加工や計算、衛星画像の素となる配列を作成
→matplotlibで可視化
→scikit-learnやpandasで機械学習用のデータセットにする　　感じです。

③コード
(1)読み込みと表示

gdal code

from osgeo import gdal #gdalはosgeoというライブラリの中にあるため、fromでimport
import numpy as np 

path="\your_path"
file_name="\data.tif"
dat=gdal.Open(path+file_name).ReadAsArray()

ReadAsArray()により、衛星画像をnumpyのarray配列に変換する。

#衛星画像の可視化
plt.figure()
plt.imshow(dat,cmap="Reds") #cmapは画像を何色で表現するか。検索すると種類が見れる
plt.colorbar()
plt.axis("off") #図の軸をすべてなくすことで衛星画像を見やすくする。
plt.show()

数行で単純な可視化が可能です。しかし、衛星画像は7バンドなのでTrue colorにしたければ工夫が必要です。

(2)衛星画像の作成の仕方
gdalとnumpyを使って、配列からラスターデータを作成する。

origin=gdal.Open(original data)#originは、空間情報を指定するためのもの
driver=gdal.GetDriverByName("GTiff") #作成するファイル形式を指定
tif=driver.Create(path+file_name+".tif",
                   origin.RasterXSize,#ラスターデータのX方向の大きさ
                   origin.RasterYSize,#ラスターデータのY方向の大きさ
                   1, #ラスターのバンド数、multi-band raster dataなら1以上を指定する。
                    gdal.GDT_UInt16 #ラスターデータのデータ型の決定)
#このデータを作成した時点で、ファイルパスのところにからのラスタデータが作成される 

#データ型は大きい値を扱うときに重要な要素である。
   
flame=np.full((XSize,YSize,1),0).T
#flameは0で満たされている配列、tifと同じshapeにする必要がある。
#あとはflameに値を以下のように渡した後で、tifのほうにflameを書き込めばraster dataの完成
#flame[band,X,Y]と指定する。mono-bandならflame[X,Y]となる
flame[:,0,0]=1
flame[:,0,3]=1#などのように
#GetRasterBand(number)で書き込むバンドを指定
#WriteArray(array)で同じshapeの配列をraster dataに書き込む
tif.GetRasterBand(n+1).WriteArray(flame)

#しかし、このままでは位置情報がないため、それを指定する。
tif.SetProjection(origin.GetProjection())#CRS情報のセット
tif.SetGeoTransform(origin.GetGeoTransform())#位置情報のセット
tif.SetMetadata(origin.GetMetadata())#メタデータのセット
tif.FlushCache()#保存
tif=None #ラスターデータを閉じることでデータの出力が終了する

(3)衛星画像の一次元化

d1=flame.flatten()#一行で完了