とらりもんHOME  Index  Search  Changes  Login

GDAL_python

How to use GDAL library on python. pythonでGDALライブラリを使うためのメモ。

Official tutorials: https://gdal.org/tutorials/raster_api_tut.html

Landsat画像(GeoTIFF)の読み込み例 / An example of reading GeoTIFF file of Landsat

from osgeo import gdal, gdal_array
ds=gdal.Open(fn_head+"B1.TIF", gdal.GA_ReadOnly)
cwidth = ds.RasterXSize               # "coarse lattice width"
cheight = ds.RasterYSize              # "coarse lattice height"
Lorig=np.empty((cheight, cwidth, 6))  # "Landsat original"

ds=gdal.Open(fn_head+"B1.TIF", gdal.GA_ReadOnly)
Lorig[:,:,1]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open(fn_head+"B2.TIF", gdal.GA_ReadOnly)
Lorig[:,:,2]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open(fn_head+"B3.TIF", gdal.GA_ReadOnly)
Lorig[:,:,3]=ds.GetRasterBand(1).ReadAsArray()

マルチバンド画像(GeoTIFF)の読み込み / How to read a multi-band GeoTIFF data

from osgeo import gdal, gdal_array
ds=gdal.Open("GPMKuPR_composite/2015_xx/GPMKuPR_2015_xx_sigma0_dB_ave.tif", gdal.GA_ReadOnly)

x=ds.ReadAsArray() # 全バンドのラスターデータを一気に読み込み / Reading all the raster data at once
x=ds.GetRasterBand(3).ReadAsArray() # 3番めのバンドのデータが得られる。 / You can get the third band.

ds.GetDescription() # ファイル名が得られる。 / You can get the filename.
ds.GetGeoTransform() # 座標系関係のパラメータが得られる。 / You can get coordinate info.
ds.GetProjection() # 投影法関係のパラメータが得られる。 / You can get projection info.

ds.GetRasterBand(3).GetDescription() # 3番めのバンドの名前が得られる。 / You can get the name of the band 3.

ポイント:

  • 画像データは, osgeo.gdal.Datasetクラスosgeo.gdal.Bandクラス(C++のGDALRasterBandクラスに相当)で扱う。
  • gdal.Openした直後のクラスはosgeo.gdal.Datasetクラス。いくつものバンドが一緒になって入っている。この時点でReadAsArray()メソッドを当てると, 全バンドのデータがnumpy.ndarrayとして得られる。
  • osgeo.gdal.Datasetクラス(のインスタンス)にGetRasterBand()メソッドを当てることで, 個々のバンドをosgeo.gdal.Bandクラス(のインスタンス)として扱うことができる。この時点でReadAsArray()メソッドを当てると, そのひとつのバンドのデータがnumpy.ndarrayとして得られる。

マルチバンド画像(GeoTIFF)の作成 ... 既にあるファイルをお手本にしてやる場合

from osgeo import gdal, gdal_array
filename='GPMKuPR_composite/2018_05/GPMKuPR_2018_05_gamma0c_dB_ave.tif'
src_ds = gdal.Open(file, gdal.GA_ReadOnly)     # お手本になるGeoTIFFファイルを読む。dsはdatasetの意味(osgeo.gdal.Datasetクラスのインスタンス)
driver = gdal.GetDriverByName("GTiff")         # GeoTIFFファイルを作るドライバー。
dst_ds = driver.CreateCopy("newfile.tif", ds)  # 新しいファイルに対応するインスタンス(osgeo.gdal.Datasetクラス)を作る。これを実行した時点でファイルができるが, データはまだNULLで埋まってる。
b = dst_ds.GetRasterBand(1)  # バンド1を呼び出す。中身はまだ空っぽ(NULL)なのだが, とにかく呼び出す。
x=b.ReadAsArray()  # バンド1のデータをnp.ndarrayに読み込む。ぜんぶNULLなのだが, 縦横の形を知りたいので。既に知ってたら不要。
x[:,:]=1            # 書き込みたいデータを作る。
b.WriteArray(x)    # バンド1のデータを, 書き込みたいデータ(この場合はx)で埋め直す。同様のことを他のバンドについても行う。bという変数名を使いまわしても大丈夫!
dst_ds.SetMetadataItem("category1", "water")   # メタデータを, 自由なキーワードで追加できる。
dst_ds.SetMetadataItem("category2", "urban")
dst_ds=None         # 新しいファイルに対応するインスタンスを破壊する。これによって, 新しいデータがファイルに書き込まれる。

既存のファイルを修正

たとえば, パラメータGeoTransformだけを変更(書き換え)したい時,

ds=gdal.Open("filename.tif", gdal.GA_Update)
ds.SetGeoTransform((-180.0, 0.05, 0.0, 40.0, 0.0, -0.05))
ds=None

とすればよい。あっさりできるよ!!

Last modified:2019/08/28 20:36:42
Keyword(s):
References:[GDAL]