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
Keyword(s):
References:[GDALとpythonでGeoTIFF画像を作成] [GDAL]