とらりもんHOME  Index  Search  Changes  Login

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

GDALライブラリを使う場合

基本はgdalライブラリでデータを読み込み、numpyの配列に変換して処理する。

参考になるページ

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)
Last modified:2018/08/21 07:27:37
Keyword(s):
References:[memo_Katagi] [とらりもんHOME]