PythonによるGeoTiff画像の読み込みと書き出し

(2019/07)一部リンクを更新

ノートの整理。
コードはインタラクティブシェルを想定。

#画像のピクセル配列データのnumpy arrayとして読み込み

read_tif.py

import numpy as np
from osgeo import gdal, gdalconst, gdal_array

src = gdal.Open('inputfilepath', gdalconst.GA_ReadOnly) # tifの読み込み (read only)
type(src) # "osgeo.gdal.Dataset"

src.RasterXSize # 水平方向ピクセル数
src.RasterYSize # 鉛直方向ピクセル数
src.RasterCount # バンド数

b1 = src.GetRasterBand(1).ReadAsArray() # 第１バンド numpy array
b2 = src.GetRasterBand(2).ReadAsArray() # 第２バンド numpy array
b3 = src.GetRasterBand(3).ReadAsArray() # 第３バンド numpy array

dtid = src.GetRasterBand(1).DataType # 型番号 (ex: 6 -> numpy.float32)
gdal_array.GDALTypeCodeToNumericTypeCode(dtid) # 型番号 -> 型名 変換

# others...

src.GetRasterBand(1).GetNoDataValue()
src.GetRasterBand(1).GetMinimum()
src.GetRasterBand(1).GetMaximum()
src.GetRasterBand(1).GetScale()
src.GetRasterBand(1).GetUnitType()

最後５項目については、自分の所持しているファイルからは何も返さなかった。
事前にメタデータのような形で番地に格納された値を参照しているのかもしれない。
もちろん、ReadAsArray()後のnumpy処理でなんとかなるケースも多い。

#基礎情報の取得

read_info.py

src.GetDescription() # ファイル名

src.GetGeoTransform() # 座標に関する６つの数字 (下記参照)
src.GetProjection() # 座標系情報

GetGeoTransform()で出力される数列の意味は、
　[始点端x座標（経度）,
　 x方向（西東）解像度,
　回転,
　始点端y座標（緯度）,
　回転,
　y方向（南北）解像度（北南方向であれば負）] 。

#画像メタデータの取得
Tiffにおけるメタデータの概念やアイテム等の検索は、例えば以下を参照。主要なメタデータについてはここに書かれているタグ名を用いてメタデータを取得する。
TIFFのフォーマット（その１） | JProgramer
CGファイル概説　第5章　第1節
GTiff -- GeoTIFF File Format

read_meta.py

src.GetMetadata() # 辞書型のメタデータ配列
src.GetMetadataItem('itemname') # 項目指定 (ex: 'TIFFTAG_XRESOLUTION', 'TIFFTAG_DATETIME', ..)

#GeoTiff画像の出力
GTiff -- GeoTIFF File Formatに従う。特に、オプションとして指定されているものも含め、しっかりと設定をしないとその後の画像使用時に正しい動作がなされないことがあるので、注意が必要。例えば、出力した画像を何かのソフトを使って処理をするような場合、そのソフトがメタデータ（タグ）のいくつかを基本情報として参照している可能性がある。その情報に欠損があると、期待した動作をしなくなる可能性があり、事前にソフトに関する理解を十分にする以外はその都度その都度原因を特定しなくてはならない。（下記コードでPhotometricInterpretationタグにRGBを入力しているのは、使用ソフト内でカラー画像で描画されると自分が期待したのに対し、第一層のバンド値を用いたモノクロ表示になる不具合があったため。）

write_data.py

from osgeo import osr # 空間参照モジュール

dtype = gdal.GDT_Float32 #others: gdal.GDT_Byte, ...
band = 3 # バンド数
output = gdal.GetDriverByName('GTiff').Create('outputfilepath', xsize, ysize, band, dtype, options = ['PHOTOMETRIC=RGB']) # 空の出力ファイル

output.SetGeoTransform((ul_x, h_res, 0, ul_y, 0, -v_res)) # 座標系指定
srs = osr.SpatialReference() # 空間参照情報
srs.ImportFromEPSG(32648) # WGS84 UTM_48nに座標系を指定
output.SetProjection(srs.ExportToWkt()) # 空間情報を結合

output.GetRasterBand(1).WriteArray(b1)   # 赤バンド書き出し（b1はnumpy 2次元配列）
output.GetRasterBand(2).WriteArray(b2)   # 緑バンド
output.GetRasterBand(3).WriteArray(b3)   # 青バンド
output.FlushCache()                     # ディスクに書き出し
output = None