📌

GeoPackageの使い方（ラスタデータ編）

2020/12/09に公開

1件

はじめに

GeoPackageでラスタファイルを扱う方法を調べてみました。

そもそもラスタってなんぞ？という方は
こちらのboiledorange73さんの記事を参照するといいと思います。

PostGIS ラスタ入門

用意しておくと便利なもの

GDAL
データ変換で使います。
こいつが肝です。
QGIS3.10
データの確認に便利です。
インストールするとGDALも入りますのでGDALのバージョンにこだわりが無ければこちらで。
spatialite-gui
GeoPackageの中身を覗くのにあると便利です。
名前の通りSpatialite用GUIツールですがGeoPackageでも利用できます。

GeoPackageにGeoTIFFデータを格納

サンプルデータを取得

サンプルデータとして使うのは「PostGISラスタ入門」と同じく農研機構のラスタデータ配信サービスの「切り取りアプリ」を使って、地図アプリケーション上から取得しましょう。

ここで使用するデータは http://aginfo.cgk.affrc.go.jp/rstprv/rstext.html.ja?lon=133.38388888889068&lat=34.385666666666815&zoom=12&layers=BTT&clon=133.3838888888889&clat=34.385666666666665&anchor=cc&iw=600&ih=400&it=image/tiff から得ています。
出典: 農研機構 (http://aginfo.cgk.affrc.go.jp/)

ダウンロードしたらdemext.tif_からdemext.tifにリネームします。

GDALを使って内容を確認しておきます。

D:\>gdalinfo demext.tif
Driver: GTiff/GeoTIFF
Files: demext.tif
Size is 600, 400
Coordinate System is:
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["JGD2000",
            DATUM["Japanese Geodetic Datum 2000",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            USAGE[
                SCOPE["unknown"],
                AREA["Japan"],
                BBOX[17.09,122.38,46.05,157.65]],
            ID["EPSG",4612]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["JGD2000 to WGS 84 (1)",
        VERSION["EPSG-Jpn"],
        METHOD["Geocentric translations (geog2D domain)",
            ID["EPSG",9603]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        USAGE[
            SCOPE["Approximation at the +/- 1m level."],
            AREA["Japan"],
            BBOX[17.09,122.38,46.05,157.65]],
        ID["EPSG",1826]]]
Data axis to CRS axis mapping: 2,1
Origin = (133.350555555559993,34.407888888888998)
Pixel Size = (0.000111111111100,-0.000111111111113)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 133.3505556,  34.4078889) (133d21' 2.00"E, 34d24'28.40"N)
Lower Left  ( 133.3505556,  34.3634444) (133d21' 2.00"E, 34d21'48.40"N)
Upper Right ( 133.4172222,  34.4078889) (133d25' 2.00"E, 34d24'28.40"N)
Lower Right ( 133.4172222,  34.3634444) (133d25' 2.00"E, 34d21'48.40"N)
Center      ( 133.3838889,  34.3856667) (133d23' 2.00"E, 34d23' 8.40"N)
Band 1 Block=600x3 Type=Float32, ColorInterp=Gray
  NoData Value=-9999

緯度経度や座標系、解像度などワールドファイルが埋め込まれたGeoTiffファイルだということが判りますね。

ちなみにGDALの環境が無い方はQGIS3系をインストールしてOSGeo4W Shellを使うと
GDALのPATHが通っていて楽ちんですよ。

demext.tifをQGISに載せるとこんな感じに表示されます。

GeoTIFFの変換方法

GeoTIFFファイルをGeoPackageに変換します。

gdalのGeoPackageラスタ関連コマンドリファレンスはこの辺り

D:\>gdal_translate -of GPKG demext.tif dem.gpkg
Input file size is 600, 400
0...10...20...30...40...50...60...70...80...90...100 - done.

指定したdem.gpkgファイルが生成されました。
さっそくspatialite_guiで覗いてみましょう。

demというテーブルにデータが格納されていることが判ります。
テーブル名を省略するとgpkgのファイル名のテーブルが作られます。

GeoTIFFの追加方法

ではdemext.tifを先ほどのgpkgファイルに別のテーブル名を指定しつつ追加してみましょう。

データの追加は
-co APPEND_SUBDATASET=YES
テーブル名の指定は
-co RASTER_TABLE=<テーブル名>

ちなみに-coはCreateOptionの略のようです。

D:\>gdal_translate -of GPKG demext.tif -co APPEND_SUBDATASET=YES -co RASTER_TABLE=demext dem.gpkg
Input file size is 600, 400
0...10...20...30...40...50...60...70...80...90...100 - done.

spatialite_guiの表示をリフレッシュすると指定したdemextテーブルが追加されています。
demextのSELECT結果を確認しても同じデータっぽい感じですね。

先ほどと同じようにGeoPackageからQGISにデータを取り込んで確認します。

正しく取り込まれていることが判りますね。

データベース構造

取り込まれたrasterデータはどのように格納されているのか見ていきます。

GeoPackageのrasterデータ構造

SELECT * FROM demext;

結果は以下のようになりました。

id	zoom_level	tile_column	tile_row	tile_data
1	2	0	0	BLOB sz=161683 TIFF image
2	2	1	0	BLOB sz=32641 TIFF image
3	2	0	1	BLOB sz=50243 TIFF image
4	2	1	1	BLOB sz=3695 TIFF image

まず最初に取り込んだGeoTIFFがタイルに変換されていることが判りますね。
そして元データの解像度から自動的に適切なズームレベルが選択されています。
画像データはバイナリデータで格納されています。
保存される画像形式はリファレンスに記載されています。
※超訳なので正しくは原文を読んでね。

デフォルトでは、GDALはPNGタイルとJPEGタイルを組み合わせて使用します。PNG形式は、入力データにアルファチャネルを含む、またはラスターの右端または下端でタイルが部分的になる場合に使用され、JPEG形式は完全に不透明なタイルの場合に使用される。

えーっと、バンド値に整数以外の値を含んでいる場合にTIFFが選択されるってのをどこかで読んだ気がするのですが、忘れてしまった…これもあとから追記します。

GeoPackageのrasterテーブル構造

次にrasterテーブルの構造を確認します。

CREATE TABLE "demext" (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    zoom_level INTEGER NOT NULL,
    tile_column INTEGER NOT NULL,
    tile_row INTEGER NOT NULL,
    tile_data BLOB NOT NULL,
    UNIQUE (zoom_level, tile_column, tile_row)
)

AUTOINCREMENTのidがPRIMARY KEY制約に指定されていますね。
ユニーク制約にズームレベルが、縦タイル番号、横タイル番号に掛ることで整合性を保っています。

GeoPackageのraster系空間演算

そもそもGeoPackageに空間演算関数はないのですが…
Spatialiteで対応していないか調べてみても残念ながら見つかりませんでした。
単なるデータストアとしてしか扱えないようです。

もしも何か知ってる人がいたら教えてもらえたら嬉しいです。

GeoPackageに背景図用rasterを格納

rasterデータには背景図としての利用方法もあります。
背景図ではピクセル毎のデータ精度はそこまで求められず、見た目が綺麗でデータ容量が軽い方が使い勝手が良いとされることが多いように思います。

サンプルデータの取得

愛知県半田市が公開しているオープンデータの航空写真を使ってみます。

図郭番号で、07ND942、07ND951、07ND944、07ND953を使います。
ワールドファイルとJPEG形式の画像をダウンロードします。

QGISに載せるとこんな感じになります。
（公開サイトに座標系などが記載されていないのでEPSG:6675で読み込んでいます。）

複数ファイルを１つのテーブルに取り込む

各タイル事に別テーブルで格納することもできますが、使い勝手が悪いので1つのテーブルにまとめて入れておくと便利です。

PostGISでいうところのraster2pgsqlの-aオプションのようなものを探したのですが見つけられませんでした。
流石にないということはないと思いますので後で調べて追記します。

仕方ないのでここではヴァーチャルラスタを作って登録する方法を試しましょう。

vrtの作成

gdalbuildvrtでヴァーチャルラスタファイルを作成します。

D:\handa>gdalbuildvrt -a_srs EPSG:6675 handa.vrt *.jpg
0...10...20...30...40...50...60...70...80...90...100 - done.

vrtファイルにはデータソースとの相対PATHが埋め込まれるので作る場所は注意してください。
上の例だと同一ディレクトリ内のファイルを参照します。

vrtの変換方法

データソースをvrtファイルに変更するだけです。

D:\handa>gdal_translate -of GPKG handa.vrt -co RASTER_TABLE=ortho handa.gpkg
Input file size is 20000, 15000
0...10...20...30...40...50...60...70...80...90...100 - done.

GeoPackageからQGISにデータを取り込んで確認します。

オプションで指定したorthoテーブルにまとめて入っていることが確認できました。

データ構造を覗いてみる

SELECT * FROM ortho;

結果は以下のようになりました。

id	zoom_level	tile_column	tile_row	tile_data
1	7	0	0	BLOB sz=667 JPEG image
2	7	1	0	BLOB sz=667 JPEG image
3	7	2	0	BLOB sz=667 JPEG image
4	7	3	0	BLOB sz=667 JPEG image
5	7	4	0	BLOB sz=667 JPEG image
6	7	5	0	BLOB sz=667 JPEG image
・
・
・
4658	7	75	58	BLOB sz=71921 PNG image
4659	7	76	58	BLOB sz=68026 PNG image
4660	7	77	58	BLOB sz=70489 PNG image
4661	7	78	58	BLOB sz=11295 PNG image

orthoテーブルのレコード数は4661です。

今回はズームレベル7、画像形式はJPEG形式とPNG形式の混在で登録されていることが判ります。
基本はJPEG、端っこはタイル全体に画像が無いので透過できるPNGになりました。

オーバービュー生成

取り込んだのは半田市の一部（たった4メッシュ）ですがQGISでレイヤを全域表示すると表示に時間がかかったと思います。
呼び出しているのは4661レコード分（トータル20000 * 15000ピクセル）の画像データになりますので重いです。
しかもQGISの画面サイズに縮小処理されてレンダリングしているのでもったいないですね。

そういうときにはオーバービューを生成することで高速化できます。

gdaladdoコマンドのリファレンス
https://gdal.org/programs/gdaladdo.html

D:\handa>gdaladdo -r cubic GPKG:handa.gpkg:ortho 2 4
0...10...20...30...40...50...60...70...80...90...100 - done.

ズームレベルにどんな変化があったのか調べてみます。

SELECT zoom_level FROM ortho GROUP BY zoom_level;

zoom_level
5
6
7

5と6のズームレベルが増えました。

オーバービューを追加

追加で増やす場合はコマンド末尾のlevelsパラメータを変えて再実行

D:\handa>gdaladdo -r cubic GPKG:handa.gpkg:ortho 8 16
0...10...20...30...40...50...60...70...80...90...100 - done.

SELECT zoom_level FROM ortho GROUP BY zoom_level;

zoom_level
3
4
5
6
7

orthoテーブルのレコード数は6261件です。

オーバービューを削除

消したい場合は-cleanオプションを付けて実行

D:\handa>gdaladdo -r cubic GPKG:handa.gpkg:ortho -clean
0...10...20...30...40...50...60...70...80...90...100 - done.

SELECT zoom_level FROM ortho GROUP BY zoom_level;

zoom_level
7

レベルを明示せずオーバービューを生成

コマンド末尾のlevels（2 4 8 16）を省略するとタイル幅が256pxよりも小さくなるまでズームレベルが生成されます。

D:\handa>gdaladdo -r cubic GPKG:handa.gpkg:ortho
0...10...20...30...40...50...60...70...80...90...100 - done.

SELECT zoom_level FROM ortho GROUP BY zoom_level;

zoom_level
0
1
2
3
4
5
6
7

orthoテーブルのレコード数は6270件です。

ズームレベルを増やしてもレコード件数や容量は微増しかしませんので普段使いはこれでいいような気がします。

QGISで確認

ではこの状態でQGISにレイヤを追加してみましょう。
初期状態とは打って変わってさくさく動いてくれたと思います。

GeoTIFFファイルのピラミッディングではなくXYZTilesのようなオーバービューが構築されました。
GeoTIFFよりも軽量で、XYZTilesの欠点であるファイル数が大量になってしまうこともないのでとても使い勝手良いですね。

Spatialite_guiで確認

カラータイルであればSpatialite_guiからもデータを確認することができます。
テーブル開いて、見たいデータの上で右クリック-[BLOB explore]

Imageタブに切り替えると

確認できました。

GeoPackageのrasterテーブル削除

rasterテーブルの削除方法

rasterテーブルはgpkg_contentsやgpkg_metadata_reference等々に関連データを持っていますので対象テーブルをばっさりDropしてしまうと残骸が残ってしまうためお勧めしません。
更にベクトルテーブル削除用のDropGeoTable関数は有効ではありませんでした。
たぶん何かあると思うのですが見つけられなかった。（きっとあるはず！）

上記のサイトでは以下のようなコマンドで削除すると書いてありました。

ogrinfo -sql "drop table raster2" droptest.gpkg
ogrinfo -sql "delete from gpkg_contents where table_name='raster2'" droptest.gpkg
ogrinfo -sql "delete from gpkg_metadata_reference where table_name='raster2'" droptest.gpkg --config OGR_GPKG_FOREIGN_KEY_CHECK OFF
ogrinfo -sql "delete from gpkg_tile_matrix where table_name='raster2'" droptest.gpkg --config  OGR_GPKG_FOREIGN_KEY_CHECK OFF
ogrinfo -sql "delete from gpkg_tile_matrix_set where table_name='raster2'" droptest.gpkg --config OGR_GPKG_FOREIGN_KEY_CHECK OFF
ogrinfo -sql "VACUUM" droptest.gpkg

今のところQGISのブラウザパネルからテーブルを削除するときれいに削除されることを確認したのでそれが手っ取り早いかと思われます。