$ 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でイケるはず...多分
$ 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
$ 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
# 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 は警告を出す。
# /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 |
# 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_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 ./
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
<!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>
$ 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));
$ /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
SELECT AddRasterConstraints('alos2_tbl','rast');
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
<!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>