とらりもんHOME  Index  Search  Changes  Login

とらりもん - Pythonで衛星データをいじる Diff

  • Added parts are displayed like this.
  • Deleted parts are displayed like this.

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 "if" "ifile" as the instance which has been opened by rasterio.open
# (既存のGeoTifファイルと同じ投影法・範囲で保存したいとする。ififileを既存の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=if.crs, transform=if.transform)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|https://rasterio.readthedocs.io/en/latest/index.html]]
**[[Quick Start|https://rasterio.readthedocs.io/en/latest/quickstart.html]]
**[[Example|https://github.com/mapbox/rasterio/tree/master/examples]]
**[[Switching from GDAL’s Python bindings|https://mapbox.github.io/rasterio/switch.html]] ... 後述するGDALライブラリとの違いを述べている。

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


!! 参考になるページ
*[[Python GDAL/OGR Cookbook|https://pcjericks.github.io/py-gdalogr-cookbook/]] ... これとgdalのAPIを眺めれば何でもできそう。

*[[Python GDAL: Write new raster using projection from old|http://gis.stackexchange.com/questions/57005/python-gdal-write-new-raster-using-projection-from-old]]
*[[Obtain Latitude and Longitude from a GeoTIFF File|http://stackoverflow.com/questions/2922532/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)