GIS PostGIS Raster画像

目次

FrontPage
・さくらVPS
・Fedora13
・SuSe10
・Docker
・Ansible
・Java
・Scala
・Python
・Ruby
・Lisp
・Computer
・GIS
・HTML
・Culture
・Random
・Link

訪問者

total: 449
today: 3
yesterday: 2
now: 9

更新

最新の10件

人気の10件

MenuBar

これは何? †

• 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 の凸なジオメトリ

コピペ用

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) に当てはめる
• 参考実装 (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: {
                'LAYERS' : 'wmap',
              },
              wrapX: true
            })
          }),
          new ol.layer.Tile({
            source: new ol.source.TileWMS({
              url : 'http://localhost:10080/cgi-bin/mapserv?map=/opt/maps/raster2.map',
              params: {
                // these are simply added to http-get parameter
                'LAYERS' : 'tokyo',
              },
              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>
  

• 実行イメージ
  • NODATA が透明色で描画されている
    
  • お台場らへん
    
• このままだと、Webなどリアルタイム性が必要な用途には耐えられないほど遅いので GIS Mapcache mapcache_seed でキャッシュを温める につづく

GIS