GIS PostGIS Raster画像

2024-10-122024
78% (285.6/366)
October
37.7% (11.6/31)
Week 41
81.2% (5.6/7)
Day 12 Sat
68.9% (16.5/24)

訪問者

total: 8845
today: 1
yesterday: 1
now: 1

↑

更新

MenuBar

これは何?
サンプルデータのダウンロード
Raster テーブルを作る
Table へ Raster 画像を登録する
ラスタ制約をつける
GDAL でラスタ画像を切り出してみる
Mapserver から Postgis の raster テーブルを参照する
衛星写真を重ね合わせる

これは何? †

PostGIS に Raster画像を入れて Mapserver から参照したい
Raster画像とは、地理空間情報がついた画像
今まで、地点や区画などの Vector データを扱ってきたけど、PostGIS は画像も取り扱うことができる

実行環境 vagrant + ansible

$ git clone https://github.com/kagyuu/ansible_gis.git
$ git clone https://github.com/kagyuu/ansible_common.git
$ cd ansible_gis
$ vagrant up
$ ansible-playbook gis.yml
$ vagrant ssh

でイケるはず...多分

↑

サンプルデータのダウンロード †

とりあえず Natural Earth の 1:50m の世界地図をダウンロード http://www.naturalearthdata.com/downloads/50m-raster-data/50m-cross-blend-hypso/

$ axel -a http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/raster/HYP_50M_SR_W.zip
$ unzip HYP_50M_SR_W.zip
$ ls -lh ./HYP_50M_SR_W
 合計 167M
 -rw-r--r-- 1 vagrant vagrant  29K 11月  8  2012 HYP_50M_SR_W.README.html
 -rw-r--r-- 1 vagrant vagrant    5 10月 18  2014 HYP_50M_SR_W.VERSION.txt
 -rw-r--r-- 1 vagrant vagrant  147  9月 19  2012 HYP_50M_SR_W.prj
 -rw-r--r-- 1 vagrant vagrant  173 12月 22  2009 HYP_50M_SR_W.tfw
 -rw-r--r-- 1 vagrant vagrant 167M 10月 18  2014 HYP_50M_SR_W.tif

tif (GeoTiff?) の属性

$ gdalinfo HYP_50M_SR_W.tif
Driver: GTiff/GeoTIFF
Files: HYP_50M_SR_W.tif
       HYP_50M_SR_W.tfw
Size is 10800, 5400
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-179.999999999999972,90.000000000000000)
Pixel Size = (0.033333333333330,-0.033333333333330)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2014:10:18 09:46:12
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CC 2014 (Macintosh)
  TIFFTAG_XRESOLUTION=342.85699
  TIFFTAG_YRESOLUTION=342.85699
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
Center      (  -0.0000000,   0.0000000) (  0d 0' 0.00"W,  0d 0' 0.00"N)
Band 1 Block=10800x1 Type=Byte, ColorInterp=Red
Band 2 Block=10800x1 Type=Byte, ColorInterp=Green
Band 3 Block=10800x1 Type=Byte, ColorInterp=Blue

画像の大きさ 10800px × 5400px / RGBカラー画像
測地系 EPSG:4326 (WGS84)
経度 -180° 〜 180°緯度 -90°〜90°
1°四方 = 30px 四方

とりあえず世界地図で練習
衛星写真 geotiff の扱いは後でやる
- nodata の画素は無視(透明色にして)欲しい
- 同じ地点で違う撮像日の画像が有るので、古い順に重ねて、新しいのが上になるように表示して欲しい

↑

Raster テーブルを作る †

# psql -U gis -h localhost -W
ユーザ gis のパスワード:
psql (9.6.2)
"help" でヘルプを表示します.

gis=> \d
                  リレーションの一覧
 スキーマ |       名前        |     型     |  所有者
----------+-------------------+------------+----------
 public   | geography_columns | ビュー     | postgres
 public   | geometry_columns  | ビュー     | postgres
 public   | raster_columns    | ビュー     | postgres
 public   | raster_overviews  | ビュー     | postgres
 public   | spatial_ref_sys   | テーブル   | postgres
 topology | layer             | テーブル   | postgres
 topology | topology          | テーブル   | postgres
 topology | topology_id_seq   | シーケンス | postgres
(8 行)

gis=> CREATE TABLE wmap_tbl (
gis(>   id bigserial primary key,
gis(>   rast raster,
gis(>   filename text
gis(> );
CREATE TABLE
gis=> CREATE INDEX idx_wmap_tbl_id ON wmap_tbl USING GiST (ST_ConvexHull(rast));
CREATE INDEX
gis=> \q

index に指定した ST_ConvexHull?(raster) は、raster の凸なジオメトリ

IMG SEARCH ERR for "ST_ConvexHull"

コピペ用

CREATE TABLE wmap_tbl (
  id bigserial primary key,
  rast raster,
  filename text
);
CREATE INDEX idx_wmap_tbl_id ON wmap_tbl USING GiST (ST_ConvexHull(rast));

※ テーブル名は小文字にする。ラスタ制約テーブルのテーブル名列に、テーブル名が小文字化されて入るので、テーブル名は小文字で入れる。実害はないかもしれないけれども、少なくとも大文字テーブル名にラスタ制約を入れても gdalinfo は警告を出す。

↑

Table へ Raster 画像を登録する †

SQL文を作成

# /usr/pgsql-9.6/bin/raster2pgsql -a -F -s 4326 -t 30x30 -Y HYP_50M_SR_W.tif wmap_tbl > raster.sql

-a	INSERT文のみ
-F	filename列にファイル名を格納
-s <srid>	空間参照系
-t <tile size>	ここで指定した 30x30 は、緯度経度 1°四方。ここで指定したサイズが一度に吐き出される最小のサイズになる。いくつでも良いけど余りが出ないようにする
-Y	Use COPY statement instead of INSERT

SQL文を実行する

# psql -U gis -d gis -h localhost -f raster.sql

↑

ラスタ制約をつける †

# gdalinfo 'PG:host=localhost dbname=gis user=gis password=password table=wmap_tbl mode=2'
Warning 1: Cannot find (valid) information about public.wmap_tbl table in raster_columns view. The raster table load would take a lot of time. Please, execute AddRasterConstraints PostGIS function to register this table as raster table in raster_columns view. This will save a lot of time.
Driver: PostGISRaster/PostGIS Raster driver
Files: none associated
Size is 10800, 5400
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.033333333333312,-0.033333333333312)
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
Center      (  -0.0000000,   0.0000000) (  0d 0' 0.00"W,  0d 0' 0.00"N)
Band 1 Block=2048x2048 Type=Byte, ColorInterp=Red
Band 2 Block=2048x2048 Type=Byte, ColorInterp=Green
Band 3 Block=2048x2048 Type=Byte, ColorInterp=Blue

"mode=2" は、テーブル全体を一枚のラスタ画像と見なす。
それはそうと、AddRasterConstraints? でラスタ制約を付けないと、検索がめちゃ遅いよという警告が出る。
ということで

# psql -U gis -h localhost -W
ユーザ gis のパスワード:
psql (9.6.2)
"help" でヘルプを表示します.

gis=> SELECT AddRasterConstraints('wmap_tbl','rast');
NOTICE:  Adding SRID constraint
NOTICE:  Adding scale-X constraint
NOTICE:  Adding scale-Y constraint
NOTICE:  Adding blocksize-X constraint
NOTICE:  Adding blocksize-Y constraint
NOTICE:  Adding alignment constraint
NOTICE:  Adding number of bands constraint
NOTICE:  Adding pixel type constraint
NOTICE:  Adding nodata value constraint
NOTICE:  Adding out-of-database constraint
NOTICE:  Adding maximum extent constraint
 addrasterconstraints
----------------------
 t
(1 行)

AddRasterConstraints?は、現在入っている raster 画像がぎりぎり満たす制約になっている。この制約外の raster 画像を追加する場合は、いったん制約を削除してから追加する

gis=> SELECT DropRasterConstraints('wmap_tbl','rast');
... raster 追加
gis=> SELECT AddRasterConstraints('wmap_tbl','rast');

↑

GDAL でラスタ画像を切り出してみる †

# gdal_translate -projwin 123 46 154 20 'PG:host=localhost dbname=gis user=gis password=password table=wmap_tbl mode=2' jp.tif
Input file size is 10800, 5400
Computed -srcwin 9090 1320 930 780 from projected window.
0...10...20...30...40...50...60...70...80...90...100 - done.

$ scp -i .vagrant/machines/gis/virtualbox/private_key -P 2222 vagrant@localhost:/home/vagrant/HYP_50M_SR_W/jp.tif ./

↑

Mapserver から Postgis の raster テーブルを参照する †

/opt/maps/raster.map

MAP
  CONFIG "MS_ERRORFILE" "/tmp/mapserver.log"
  IMAGETYPE PNG
  EXTENT -180 -90 180 90 # <Lower Left X> <Lower Left Y> <Upper Right X> <Upper Right Y>
  SIZE 256 256
  PROJECTION
    "init=epsg:4326"
  END

  LAYER
    NAME wmap
    METADATA
      "wms_title" "wmap_50m"
      "wms_srs" "EPSG:4326 EPSG:3857"
      "wms_enable_request" "*"
    END
    STATUS DEFAULT
    TYPE RASTER
    DATA "PG:host=localhost dbname=gis user=gis password=password table=wmap_tbl mode=2"
    PROCESSING "NODATA=0"
    # PROCESSING "SCALE=AUTO"
  END
END # End of MAP

参考実装 (Openlayers3)

<!doctype html>
<html lang="en">
  <head>
    <link rel="stylesheet" href="https://openlayers.org/en/v4.0.1/css/ol.css" type="text/css">
    <style>
      .map {
        height : 480px;
        width  : 800px;
        border : 1px solid black;
      }
    </style>
    <script src="https://openlayers.org/en/v4.0.1/build/ol.js" type="text/javascript"></script>
    <title>OpenLayers 3 example</title>
  </head>
  <body>
    <div id="map" class="map"></div>
    <script type="text/javascript">
      var layers = [
        new ol.layer.Tile({
          source: new ol.source.TileWMS({
            url : 'http://localhost:10080/cgi-bin/mapserv?map=/opt/maps/raster.map',
            params: {
              // these are simply added to http-get parameter
              'LAYERS' : 'wmap',
            },
            wrapX: true
          })
        })        
      ];
 
      var map = new ol.Map({
        target: 'map',
        layers : layers,
        view: new ol.View({
          center: ol.proj.transform([135.0, 35.0], 'EPSG:4326', 'EPSG:3857'),
          zoom: 4
        })
      });
    </script>
  </body>
</html>

実行イメージ

↑

衛星写真を重ね合わせる †

JAXA から ALOS-2 のサンプル画像をダウンロード (東京解像度3m)

$ axel -a ftp://ftp.eorc.jaxa.jp/pub/ALOS-2/1501sample/020_tokyo/0000008655_001001_ALOS2014410740-140829.zip
$ gdalwarp -t_srs EPSG:4326 IMG-HH-ALOS2014410740-140829-UBSL3.1GUA.tif alos2.tif
Creating output file that is 32659P x 29912L.
Processing input file IMG-HH-ALOS2014410740-140829-UBSL3.1GUA.tif.
0...10...20...30...40...50...60...70...80...90...100 - done.

今回は変な測地系 ( https://epsg.io/4338 ) の GeoTiff? なので 4326 (WGS84=緯度経度) に変換する

テーブル作成

CREATE TABLE alos2_tbl (
  id bigserial primary key,
  rast raster,
  filename text
);
CREATE INDEX idx_alos2_tbl_id ON alos2_tbl USING GiST (ST_ConvexHull(rast));

raster画像の登録

$ /usr/pgsql-9.6/bin/raster2pgsql -a -F -s 4326 -t 256x256 -Y alos2.tif alos2_tbl > alos2.sql
$ psql -U gis -d gis -h localhost -f alos2.sql

ラスタ制約忘れずに(さすがに1.8GBのデータなんで、ラスタ制約付けないとマジ遅くなる)
```
SELECT AddRasterConstraints('alos2_tbl','rast');
```

/opt/maps/raster2.map

MAP
  CONFIG "MS_ERRORFILE" "/tmp/mapserver.log"
  IMAGETYPE PNG
  EXTENT -180 -90 180 90 # <Lower Left X> <Lower Left Y> <Upper Right X> <Upper Right Y>
  SIZE 256 256
  PROJECTION
    "init=epsg:4326"
  END

  LAYER
    NAME tokyo
    METADATA
      "wms_title" "tokyo_3m"
      "wms_srs" "EPSG:4326 EPSG:3857"
      "wms_enable_request" "*"
    END
    STATUS DEFAULT
    TYPE RASTER
    DATA "PG:host=localhost dbname=gis user=gis password=password table=alos2_tbl mode=2"
    #PROCESSING "BANDS=1"
    PROCESSING "NODATA=0"
    PROCESSING "SCALE=AUTO"
  END # End of LAYER
END # End of MAP

ミソは NODATA と SCALE
NODATA=0 : 値が0の場所はデータがないものと見なす → 透明色で描画する
SCALE=AUTO : 1band 16bit データを適当に丸めて RGB(0-255,0-255,0-255) に当てはめる