Python3×地理空間データ GDAL Python API 【未完】
この記事は以下のような方を対象としています.
- Pythonで地図データ,位置データを操ってみたい
- プログラミングの基本的なことはわかるはずだ
- 地図・GISに関する知識がある
- 知識がなくても自分で調べられる
- tarokと一緒に勉強したい(笑)
- 未熟なtarokにアドバイスしたい (ありがたい)
はじめに
「Python×地理空間データ」シリーズの第3弾です.
前回ご紹介したGDALのPythonAPIについて,より詳しくまとめます.
今回,主に参考にさせていただいたWebページは以下の通りです.
なお,ここではサンプルデータを用い,可能な限り具体的にまとめていきたいと思います.
サンプルデータ
今回使用するラスタデータは,岡山県南部の数値標高モデル(DEM:Digital Elevation Model)です.
QGISで可視化すると,このような感じです.
DEMデータソース | 国土地理院 基盤地図情報 数値標高モデル10mメッシュ |
DEM凡例 |
なお,コードの実行環境は下表の通りです.
OS | Windows |
プログラミング言語 | Python 3.7.3 (in Anaconda) |
GDAL/OGR | gdal 2.3.3 |
GDALの概要
ライブラリ構造
GDALはOGRと統合され,GDAL/OGRという単一ライブラリとして提供されていると前回書きました.
Python(Anaconda)においては「conda install gdal」でそれがインストールできまして,中身は以下のような構造をしています.
インストールするときはgdalなのに,パッケージ名はosgeoなんですね.ちょっとややこい.
まあそれは置いといて,今回取り上げるモジュールは主にgdalで,少しgdal_arrayやgdalconstも登場します.
他のモジュールは別の機会で扱います.
データモデル
前回の画像を再掲します.
ラスタデータは,GDALではこの図のようなオブジェクトとして扱われるのでした.
データセットの中にラスターバンドが1つ以上含まれる構造です.
バンドはピクセル値の行列であり,これがラスタデータの主要部とも言えるでしょう.
ところで,このデータモデルでは,あたかもBandがDatasetを継承しているように錯覚してしまいます.(僕だけかもしれませんw)
しかし,実際は下図のように,両者がMajorObjectを継承しているようです.
MajorObjectは「メタデータ」を扱うための汎用的なクラスであり,DatasetやBand以外にもMajorObjectを継承しているクラスがあります.(ここでは割愛します.)
つまり,メタデータはMajorObjectで扱い,より具体的なクラスの処理は各派生クラスで実装するというデザインなんでしょう.
継承関係とかややこしい話になってしまいましたが,とりあえず先程のデータモデルを念頭においてプログラミングしていけば十分だと思います.
ラスタデータの読み込み
from osgeo import gdal
dataset = gdal.Open("hogehoge.tif")
gdal.Open()を使って読み込みます.
読込モード
gdal.Open()では,読込オンリーモード(GA_ReadOnly)と更新モード(GA_Update)を指定できます.
上記のように何も指定しなければ自動的に読込オンリーモードになります.
データセット(ファイル)の値を直接変更したい場合は,以下のように更新モードにしてください.
dataset = gdal.Open("hogehoge.tif", gdal.GA_Update)
また,更新モードで中身を変更した後,以下のようにファイルへの書き込みをしないと反映されません.
dataset.FlushCache() #方法1:書き込むだけ
dataset = None #方法2:書き込んでファイルを閉じる
読むだけで十分な場合や,既存のデータセットの情報を用いて新しいデータセットを作成する場合は,読込モードで十分です.
Dataset
help()でDatasetの中身を覗いてみたらわかるのですが,かなり多くのメソッドやプロパティがあります.
私は面倒くさがりゆえ,データモデルに関連するものを中心にまとめます.
プロパティ
各項目はアコーディオンになっているので,参照したい項目をクリックしてください.
概要
Datasetが含むバンドの数
出力
int
サンプル使用例
dataset.RasterCount
1
標高値の1バンドしか入ってません.
概要
幅(東西方向のピクセル数)
出力
int
サンプル使用例
dataset.RasterXSize
3032
概要
高さ(南北方向のピクセル数)
出力
int
サンプル使用例
dataset.RasterXSize
4063
XSize(3032)×YSize(4063)=12,319,016個のセル数ってことですね.
ビックデータだあぁ
メソッド
概要
空間参照系に関する情報を取得する.
出力
string(以下のフォーマットのいずれか)
- Well Known Text (WKT)
- Human-Readable OGC WKT
- Proj4
- OGC WKT
- JSON
- GML
- ESRI WKT
- .PRJ File
- USGS
- MapServer Mapfile | Python
- Mapnik XML | Python
- GeoServer
- PostGIS spatial_ref_sys INSERT statement
- Proj4js format
フォーマットについては「spatialreference.org」で登録・検索・閲覧ができます.
サンプル使用例
dataset.GetProjection()
'LOCAL_CS["JGD2000 / Japan Plane Rectangular CS V",
GEOGCS["JGD2000",
DATUM["unknown", SPHEROID["unretrievable - using WGS84", 6378137,298.257223563], TOWGS84[0,0,0,0,0,0,0]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
AUTHORITY["EPSG","2447"],
UNIT["metre",1, AUTHORITY["EPSG","9001"]]]'
上の結果は,入れ子構造が見やすいように少々加工してます.
ちなみに,このようなデータをWKT(well-known text)と呼びます.
WKTについては,このサイトなどに詳しく載っています.
入れ子構造の上位から見ていきましょう.
- LOCAL_CS
ローカル測地系(Local Coordinate System)であることを示しています.
地図作成や測量などにおいて共通的に使用できる経緯度や標高を定義するための基準を測地系といいます.
この測地系には,特定地域のみで適用できるローカル測地系と,世界的範囲で適用できる世界測地系があります.
サンプルデータは日本測地系であり,ローカル測地系のひとつです.
この話は,さすがのEsriさんがわかりやすいので詳細はお譲りします.
空間参照系全般については,ここがわかりやすいです.
- JGD2000 / Japan Plane Rectangular CS V
「JGD2000」は日本測地系を,「Japan Plane Rectangular CS V」は平面直角座標系の第5系であることを意味しています.
詳細は,ここやそこが詳しいです.
- GEOGCS
地理座標であることを示しています.
投影座標であれば「PROJCS」,地心から見た座標なら「GEOCCS」になります.
ここに載っています.
- DATUM
水平方向の測地基準を示します.
- SPHEROID
選択されている回転楕円体に関する情報を示します.
ここらへんの話です.
- TOWGS84
WGS84以外の測地系からWGS84に変換するためのパラメータです.
ここに詳しく載っていました.
- PRIMEM
基準にする経線を定義しています.
このデータでは「グリニッジ子午線(Greenwich)」です.
また,数値が正の値ならグリニッジより東に,負の値なら西に基準があることを示します.
ここを参考にしました.
- AUTHORITY
このサイトによるとオプション情報を格納できる場所とのことですが,ほとんどの場合はEPSGコードが入っています.
EPSGコードとは,各座標参照系に割り振られたIDです.
このサンプルの場合は「2447」 が「JGD2000 / Japan Plane Rectangular CS V」に対応しています.
よく使う空間参照系のEPSGは覚えておくと便利です.
日本だけならここが詳しいです.
EPSGそのものについてはこのサイトに載ってます.
- UNIT
単位です.
「metre」ならばメートル単位の距離,「degree」なら角度です.
概要
空間参照系に関する情報を設定します.
入力
String(以下のフォーマットのいずれか)
フォーマットの記述法の詳細はGetProjection()やspatialreference.orgを見てください.
サンプル使用例
# 更新モードで読み込んでいるものとする
wkt = 'GEOGCS["JGD2000",DATUM["unknown",SPHEROID["unnamed",6378137,298.2572221010042],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4612"]]'
dataset.SetProjection(wkt)
dataset.FlushCache() #変更の保存
概要
地図画像の情報(座標値など)を取得します.
出力
tuple
(num, num, num, num, num, num)
中身は以下の順番です.
- 画像左上セルのX座標
- 単一セルのXサイズ(東ー西)
- 回転度(上が北なら0)
- 画像左上セルのY座標
- 回転度(上が北なら0)
- 単一セルのYサイズ(北ー南)
画像は基本的に左上端を基準にします.
使用例
dataset.GetGeoTransform()
(-65112.075071631014,
11.418434105844657,
0.0,
-129197.91163072696,
0.0,
-11.418434105844657)
概要
地図画像の情報(座標値など)を設定します.
入力
tuple
(num, num, num, num, num, num)
中身の詳細はGetGeoTransform()を見てください.
使用例
# 更新モード
gt = (-65112.075071631014, 11.418434105844657, 0.0, -129197.91163072696, 0.0, -11.418434105844657)
dataset.SetGeoTransform(gt)
dataset.FlushCache()
概要
データセットのファイル名(フルパス)を取得します.
出力
string
使用例
dataset.GetDescription()
'C:\\hogehoge\\raster.tif'
概要
メタデータ(補助的な情報)を設定します.
入力
dict
詳細はGetMetaData()をご覧ください.
サンプル使用例
# 更新モード
dataset.SetMetadata({"TIFFTAG_DOCUMENTNAME": "metameta"})
dataset.FlushCache()
概要
データセットに含まれるバンドを取得します.
入力
int
取得したいバンドの番号(1~)です.
出力
bandオブジェクト
サンプル使用例
band1 = dataset.GetRasterBand(1)
type(band1)
osgeo.gdal.Band
概要
全てのバンドの配列部分をnumpy配列(ndarray)として取得します.
マルチバンドデータの場合,すべての配列を一度に取得すると非常に重くなるので,通常は「GetRasterBand()」で取得したバンドに対して「ReadAsArray()」を適用します.
出力
ndarray
サンプル使用例
dataset.ReadAsArray()
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]], dtype=float32)
単バンドデータなので,1つの配列が取得できました.
概要
全てのバンドのラスタデータをバイト列で取得します.
出力
サンプル使用例
Band
band = dataset.GetRasterBand(1)
Bandの中身も,一部のみのご紹介とさせていただきます.
プロパティ
メソッド
- 概 要:
- 返り値:
- サンプル使用例
gdal関連モジュール
gdal_array
bandのnumpy配列(ndarray)とラスタ―データ(Dataset)間のインターフェイスに特化したモジュールです.
from osgeo import gdal_array
# ラスタデータから配列を取得
array = gdal_array.LoadFile("hogehoge.tif")
# 配列からラスタデータを生成
ds = gdal_array.OpenArray(array)
gdalconst
gdalconstは,gdal内部で使われている専門用語とそれらに付けられた定数値(ID)の扱いに特化したモジュールです.