Handle satellite data by Python: Pythonで衛星データをいじる
履歴
2016/07/25 Jin Katagi 作成 2024/12/14 Taiga Sasagawa 最終更新
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)])
gdal.Open().ReadAsArrayとrasterio.open().read()の速度比較
対象ファイル: HLS.S30.T55TFL.2024306T010801.v2.0.B03.tif (約26 MB)
- gdalを使った場合:
In [1]: from osgeo import gdal In [2]: import time In [3]: def load_geotiff(path): ...: start = time.time() ...: data = gdal.Open(path).ReadAsArray() ...: end = time.time() ...: print(end - start) In [4]: load_geotiff("./HLS.S30.T55TFL.2024306T010801.v2.0.B03.tif") 0.6520583629608154
- rasterioを使った場合:
In [1]: import rasterio as rio In [2]: import time In [3]: def load_geotiff(path): ...: start = time.time() ...: data = rio.open(path).read(1) ...: end = time.time() ...: print(end - start) In [4]: load_geotiff("./HLS.S30.T55TFL.2024306T010801.v2.0.B03.tif") 0.21122956275939941
- 結論: rasterioを使ったほうが早い!
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]