最近衛星画像の処理をする機会が増えましたが，データ容量が大きいため，GISソフトウェアなんか使っていたらいつまでも終わらないので，GDALやpythonで処理をします。
結構煩雑なので備忘録的にまとめておきます。
ちなみにフリーで手に入るものは，
ALOS AVNIR-2 https://www.eorc.jaxa.jp/ALOS/alos-ori/index.html
LANDSAT8(他にもいろんなダウンロードサイトがあります) https://landbrowser.airc.aist.go.jp/landbrowser/index.html
等があります。
回転したgeotifのRGB合成
衛星画像は緯度に対して回転した（傾いた）geotiffがほとんどですが，gdalが対応してません。なので，worldfileを使って以下のようなに処理が必要。
rem worldfile(tfw file)を作成 gdal_translate -co "TFW=YES" band.tif out.tif rem worldfileを設定（作業終了後は必ずunsetする） set GDAL_GEOREF_SOURCES=out.tfw rem vrtファイルを作成 gdalbuildvrt -separate out.vrt bandred.tif bandgreen.tif bandblue.tif rem RGB合成したtifを作成 gdal_translate out.vrt out.tif
参考サイト
https://gis.stackexchange.com/questions/271995/how-to-get-gdaltranslate-to-create-world-file
https://stackoverflow.com/questions/52334933/merging-multiple-bands-together-through-gdal-correctly
pansharpen
LANDSAT等は，パンシャープンを行う。
gdal_pansharpen panchro.tif bandred.tif bandgreen.tif bandblue.tif out.tif
参考サイト
https://gdal.org/programs/gdal_pansharpen.html
大規模tifのmerge
基本的には使わない。gdalbuildvrtでvrtを作る。
圧縮とかbigtiffとかを指定。
gdal_merge -o out.tif -ot Uint16 -co COMPRESS=LZW -co BIGTIFF=YES file1.tif file2.tif ......
参考サイト
https://qiita.com/naru-T/items/b87c83fad05634919532
ラスタ演算
pythonがかなり楽です。
import xarray as xr # この一行でgeotiffをxarrayで読み込む d = xr.open_rasterio('*.tif') # バンドを取得，シングルでもマルチでも同じ band = d.values # こんな感じで計算 ind = (band[2] - band[4]) / (band[2] + band[4])
ちなみに結果をgeotiffで出力する場合，rioxarrayを使う。
import rioxarray xrds = xr.Dataset({'value':(['y','x'], ind)} , coords={'x':d['x'], 'y':d['y']} , attrs= d.attrs ) # うまくepsgが付かないことがある。 #xrds.rio.set_crs("epsg:****") # 2020/7/7追記 xrds = xrds.rio.write_crs('EPSG:****',inplace = True) # 投影を変換する場合 # xrds = xrds.rio.reproject("epsg:****") xrds['value'].rio.to_raster('out.tif')
tellus(https://www.tellusxdp.com/ja/)がローカルのjupyterでも使えるようになったらいいのな。
2020/5/2追記
以下のサイトにも詳しく書いてました。
