とらりもん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

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ファイルの新規作成 / How to make a GeoTiff file (1 band, latlon projection)

import gdal, ogr, os, osr
array=np.array(.....).astype(np.uint8)  # 画像データ (You create!)
WEST=135.2   # longitude of the west edge (You decide!)
NORTH=36.5   # latitude of the north edge (You decide!) 
dx=0.01      # pixel size (You decide!)
dy=-0.01     # pixel size (You decide!)
cols, rows = array.shape[1], array.shape[0]
originx , originy = WEST, NORTH
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create("output.tif", cols, rows, 1, gdal.GDT_Byte)
outRaster.SetGeoTransform((originx, dx, 0, originy, 0, dy))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(array)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(4326)
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()

マルチバンド画像(GeoTIFF)の新規作成

cols, rows, nbands = 2000, 1000, 6
SRS = osr.SpatialReference()
SRS.ImportFromEPSG(4326)
originx , originy = 140.25, 35.55
dx=0.01
dy=-0.01
outRaster = gdal.GetDriverByName('GTiff').Create("out.tif", cols, rows, nbands, gdal.GDT_UInt16)
outRaster.SetGeoTransform((originx, dx, 0, originy, 0, dy))
outRaster.SetProjection(SRS.ExportToWkt())
for i in range(nbands):
  outRaster.GetRasterBand(i+1).WriteArray(data[:,:,i])
outRaster.FlushCache()
outRaster=None

マルチバンド画像(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

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

Inserting meta data (ancillary data)

band = dst.GetRasterBand(1)
band.WriteArray(raster) # rasterは2次元の配列
band.SetMetadata(meta) # metaはdictionary. 例えば{"Offset" : 0, "Slope": "1e-5"}など
band.SetNoDataValue(65535)

取り出す時は、

GetRasterBand(1).GetMetadata()

reprojection, resampling ... gdal.Warp

dsnew=gdal.Warp("new.tif", "source.tif", srcSRS="EPSG:4612", dstSRS="EPSG:4612", 
         outputBounds=(xmin, ymin, xmax, ymax), xRes=xres, yRes=yres, resampleAlg="average")

You can find how to give option parameters here https://gdal.org/python/osgeo.gdal-module.html#WarpOptions

Last modified:2022/08/26 16:57:03
Keyword(s):
References:[GDALとpythonでGeoTIFF画像を作成] [GDAL]