GDALとpythonでGeoTIFF画像を作成
GDAL_pythonに書かれている手法を試してみたものの、目的の画像が上手く作成できなかったので、代わりに使った方法です。ipython3でやりました。(嶌田)
output.SetGeoTransformの部分を自動でできるように追記しました.(笹川)
#NDVI、赤、近赤外のバンドを持つGeoTIFF画像の作成が目標です。 #必要なパッケージのインストール #rasterioを予めインストールしておいてください import rasterio as rio import numpy as np from osgeo import gdal, gdalconst, gdal_array from osgeo import osr #各種情報の取得 with rio.open('filepath') as filename : filename.bounds filename.meta #出力例{'driver': 'GTiff','dtype': 'uint16','nodata': None,'width': 6037,'height': 4800,'count': 9, 'crs': CRS({'init': 'epsg:4326'}),'transform': (-63.85066604614258, 0.002166254445910454, 0.0, -10.0, 0.0, -0.0020833334419876337), 'affine': Affine(0.002166254445910454, 0.0, -63.85066604614258, 0.0, -0.0020833334419876337, -10.0)} #widthとheightはそれぞれ横と縦方向の解像度 #transform以下の情報は(始点端経度,西東解像度,回転,始点端緯度,回転,南北解像度) ds=gdal.Open("filepath", gdal.GA_ReadOnly) Red=ds.GetRasterBand(number_of_Redband).ReadAsArray() #number_of_Redbandには使用する画像で赤色に対応するバンドの番号 NIR=ds.GetRasterBand(number_of_NIRband).ReadAsArray() #number_of_NIRbandには使用する画像で近赤外に対応するバンドの番号 NDVI=(NIR-Red)/(NIR+Red+1) #"+1"は0で割ることの防止で入れた(あんまりよくないかも)。 #自分がやった時には、-1から1の値をとるはずのNDVIが1を大きく超えてしまったので(原因は書いた時点では不明)、以下の様な感じで無理やり1を越えてる値をNaNにしました。 原因をご存知の方がいらっしゃったら追記をお願いします‥ for i in range(ds.RasterYSize): for j in range(ds.RasterXSize): if NDVI[i,j]>1: NDVI[i,j]=np.nan else: continue #ファイルの書き出し xsize=ds.RasterXSize #水平方向ピクセル数 ysize=ds.RasterYSize #垂直方向ピクセル数 band=3 #NDVIと赤と近赤外で画像を作るので3つバンドを用意 dtype = gdal.GDT_Float64 output = gdal.GetDriverByName('GTiff').Create('output_filepath', xsize, ysize, band, dtype) ul_x=-63.85066604614258 #始点端経度 h_res=0.002166254445910454 #西東解像度 ul_y=-10.0 #始点端緯度 v_res=-0.0020833334419876337 #南北解像度 output.SetGeoTransform((ul_x, h_res, 0, ul_y, 0, v_res)) #自動で取得させたい場合にはoutput.SetGeoTransform(ds.GetGeoTransform())とすれば良い.上のul_xやh_resを打ち込まなくでもできる. srs = osr.SpatialReference() #空間参照情報 srs.ImportFromEPSG(4326) #最初に取得したEPSGコードを入れる output.SetProjection(srs.ExportToWkt()) output.GetRasterBand(1).WriteArray(NDVI) #出力画像のバンド1にNDVIの計算結果を入れる output.GetRasterBand(2).WriteArray(Red) #バンド2に赤 output.GetRasterBand(3).WriteArray(NIR) #バンド3に近赤外 output.FlushCache() #書き出し output = None
Keyword(s):
References:[WGS84とUTM]