Handle satellite data by Python: Pythonで衛星データをいじる
2016/07/25 Jin Katagi 作成
Use Rasterio (the easiest!): Rasterioを使う場合 (これがいちばん楽!)
Install Raterio
sudo apt-get install python-rasterio python3-rasterio or sudo pip3 install rasterio
Read GeoTif file: Geotifファイルの読みこみ
if=rasterio.open('filename.tif') data=if.read(1) # 1 is the band ID, which can be 2, 3, 4, ... for other bands. # note that this ID starts from 1, not from 0.
Write GeoTif file: GeoTifファイルの書き出し
# Suppose you want to make same projection and region as an existing GeoTif File with "ifile" as the instance which has been opened by rasterio.open # (既存のGeoTifファイルと同じ投影法・範囲で保存したいとする。ifileを既存のGeoTifファイルをrasterio.openで開いたオブジェクトとする) # outdataというnumpy配列に書き出したいデータが入っているとする。 of=rasterio.open('filename.tif', 'w', driver='GTiff', height=data.shape[0], width=data.shape[1], count=1, dtype=data.dtype, crs=ifile.crs, transform=ifile.transform) of.write(data, 1) of.close() # 注意!! numpy配列のshapeは, 第1パラメータは縦(line数), 第2パラメータが横(pixel数)。これ, 逆じゃね?と勘違いしそう。 # 注意!! dtypeとして, Int8は使えないことがある(gdalのせいだっけ? GeoTifのせいだっけ?)。もしint8でダメならInt16でやるべし。 data=data.astype(int16) # マルチバンドデータを書き出すとき(3バンドの例) データをnp.ndarrayの(3, line数, pixel数)の形の配列にする。それをdataという名前とする。 rasterio.openでcount=3に変更 of.write(data)
- Meanings of keywords of rasterio.open():
driver: the name of the desired format driver width: the number of columns of the dataset height: the number of rows of the dataset count: a count of the dataset bands dtype: the data type of the dataset crs: a coordinate reference system identifier or description transform: an affine transformation matrix, and nodata: a “nodata” value
- Rasterio
- Quick Start
- Example
- Switching from GDAL’s Python bindings ... 後述するGDALライブラリとの違いを述べている。
GDALライブラリを使う場合
基本はgdalライブラリでデータを読み込み、numpyの配列に変換して処理する。
参考になるページ
- Python GDAL/OGR Cookbook ... これとgdalのAPIを眺めれば何でもできそう。
- Python GDAL: Write new raster using projection from old
- Obtain Latitude and Longitude from a GeoTIFF File
tifの読み込み
画像は(height, width, band)の順に読み込まれる。
import numpy as np from osgeo import gdal from osgeo import gdal_array
特定のバンドのみ読み込む
dataset1 = gdal.Open("test.tif", gdal.GA_ReadOnly) # numpyの配列に変換 band1 = np.array(dataset1.GetRasterBand(1).ReadAsArray())
すべてのバンドを読み込む
dataset1 = gdal.Open("test.tif", gdal.GA_ReadOnly) bands = dataset1.RasterCount rows = dataset1.RasterYSize cols = dataset1.RasterXSize # numpyの配列に変換 band_all = np.array([dataset1.GetRasterBand(band + 1).ReadAsArray() for band in range(bands)])
RGB合成
tifファイル一つに対して一つのバンドのみ入っている場合、numpy配列で読み込むと二次元配列になってしまう。(ex: (x行、y列))
その場合、以下のように新たな軸をnp.newaxisで指定して、np.c_で合成する。
import numpy as np
band_rgb = np.c_[band_red[:,:,np.newaxis], band_green[:,:,np.newaxis], band_blue[:,:,np.newaxis]] band_rgb.shape # => (x,y,3)
Keyword(s):
References:[memo_Katagi] [とらりもんHOME]