とらりもんHOME  Index  Search  Changes  Login

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

   

Last modified:2019/10/30 15:38:07
Keyword(s):
References:[WGS84とUTM]