2013年11月25日月曜日

23 - OpenLayers で地図を重ねる 3 - 数値標高モデルの画像の地図を重ねる2

23-3 あきる野市の地図に数値標高モデルの画像の地図を重ねる
tokyo10m.png を EPSG:2451 の投影法に変換します。

user@debian7-vmw:~/mapdata/tokyo$ gdal_translate -of PNG -a_srs "EPSG:2451" tokyo10m.png tokyo10m-epsg2451.png
Input file size is 2310, 1467
0...10...20...30...40...50...60...70...80...90...100 - done.
user@debian7-vmw:~/mapdata/tokyo$ gdalinfo tokyo10m-epsg2451.png
Driver: PNG/Portable Network Graphics
Files: tokyo10m-epsg2451.png
       tokyo10m-epsg2451.png.aux.xml
Size is 2310, 1467
Coordinate System is:
PROJCS["JGD2000 / Japan Plane Rectangular CS IX",
    GEOGCS["JGD2000",
        DATUM["Japanese_Geodetic_Datum_2000",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6612"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4612"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",36],
    PARAMETER["central_meridian",139.8333333333333],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",NORTH],
    AXIS["Y",EAST],
    AUTHORITY["EPSG","2451"]]
Origin = (138.873693408405359,36.000160021331652)
Pixel Size = (0.000487353665300,-0.000397437712900)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 138.8736934,  36.0001600) (139d50' 5.55"E, 36d 0' 1.17"N)
Lower Left  ( 138.8736934,  35.4171189) (139d50' 5.55"E, 36d 0' 1.15"N)
Upper Right ( 139.9994804,  36.0001600) (139d50' 5.59"E, 36d 0' 1.17"N)
Lower Right ( 139.9994804,  35.4171189) (139d50' 5.59"E, 36d 0' 1.15"N)
Center      ( 139.4365869,  35.7086395) (139d50' 5.57"E, 36d 0' 1.16"N)
Band 1 Block=2310x1 Type=Byte, ColorInterp=Red
Band 2 Block=2310x1 Type=Byte, ColorInterp=Green
Band 3 Block=2310x1 Type=Byte, ColorInterp=Blue

tokyo10m-epsg2451.png の四角の座標が度数（degree）なのでこのままでは表示できません。
EPSG:2451 の四角の座標を Proj4js で計算します。

OpenLayers and Proj4js
http://trac.openlayers.org/wiki/Documentation/Dev/proj4js

ここを参考にして、Proj4js を設定します。

1 ダウンロード
a Proj4js サイト
http://trac.osgeo.org/proj4js/

の「Download」をクリックします。
b Download-Proj4js
http://trac.osgeo.org/proj4js/wiki/Download

の「Proj4j 1.1.0」の「 http://download.osgeo.org/proj4js/proj4js-1.1.0.zip」をクリックします。

2 インストール
a ダウンロードした  proj4js-1.1.0.zip を解凍します。
user@debian7-vmw:~/ダウンロード$ unzip proj4js-1.1.0.zip
b proj4js フォルダを openlayersTokyoproj/lib/ に移動します。
user@debian7-vmw:~/ダウンロード$ mv proj4js ../mapsite/openlayersTokyoproj/OpenLayers-2.13.1/lib/
c 解凍されたフォルダの proj4js/lib/def/ に定義ファイルがあるので、EPSG:4326 と EPSG:2451 がないときは作成します。
EPSG4302.js を参考にします。
データは、/usr/share/proj のものを使用します。

user@debian7-vmw:/usr/share/proj$ grep -A1 4326 epsg
<4326> +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs  <>
# Anguilla 1957
user@debian7-vmw:/usr/share/proj$ grep -A1 2451 epsg
<2451> +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs  <>
# JGD2000 / Japan Plane Rectangular CS X
--
（EPSG:4326 の定義は proj4js/lib/proj4js ありました。
'EPSG:4326': "+title=long/lat:WGS84 +proj=longlat +a=6378137.0 +b=6356752.31424518 +ellps=WGS84 +datum=WGS84 +units=degrees",）

EPSG:2451 の定義ファイル EPSG2451.js を 次の内容で proj4js/lib/defs/ に作成します。
Proj4js.defs["EPSG:2451"]= "+title=JGD2000/Japan Plane Rectangular CS IX +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

proj4js/index.html ファイルに次のように追加します。

---
<script src="lib/proj4js.js"></script>

<script src="lib/defs/EPSG27200.js"></script>
<script src="lib/defs/EPSG4272.js"></script>
<script src="lib/defs/EPSG54009.js"></script>
<script src="lib/defs/EPSG42304.js"></script>
<script src="lib/defs/EPSG25833.js"></script>
<script src="lib/defs/EPSG27563.js"></script>
<script src="lib/defs/EPSG4139.js"></script>
<script src="lib/defs/EPSG4302.js"></script>
<script src="lib/defs/EPSG31285.js"></script>
<script src="lib/defs/EPSG900913.js"></script>
<script src="lib/defs/EPSG2451.js"></script> <!-- 追加 -->

---

Webブラウザに

http://localhost/mapsite/openlayersTokyoproj/OpenLayers-2-13.1/lib/proj4js/index.html

と入力して Enter キーを押します。

「source」を「EPSG:4326 - long/lat:WGS84」、
「dest」を「EPSG:2451 - JGD2000/Japan Plane Rectangular CS IX」
にして、gdalinfo tokyo10m-epsg4326.png の出力結果の「Corner Coordinates」から、出力します。

           「source」                  「dest」
Upper Left 138.8736934, 36.0001600 -> -86517.10966018954, 443.645829377676
Lower Left 138.8736934, 35.4171189 -> -87149.5276955563, -64243.31846699882
Upper Right139.9994804, 36.0001600 -> 14978.910498662082, 30.51726815643259
Lower Right139.9994804, 35.4171189 -> 15088.38925912324, -64653.61709593148


gdalinfo tokyo10m-epsg4326.png の出力結果の「Size is 2310, 1467」から、次のコマンドで座標を変換します。
座標は左上から時計回りに記述します。

gdal_translate [-of format] [-a_srs srs_def] [-gcp pixel line easting northing] [-gcp pixel line easting northing] [-gcp pixel line easting northing] [-gcp pixel line easting northing] sourcefile outpulfile

*2014.2.23
次のコマンドでOKでした。
gdal_translate [-of format] [-a_srs srs_def] [-a_ullr ulx uly lrx lry] sourcefile outpulfile
（ul: upper left, lr: lower right）
[-gcp pixel line easting northing] は、その画像の点（ピクセルで表した位置）の座標を設定するコマンドみたいです。

user@debian7-vmw:~/mapdata/tokyo$ gdal_translate -of PNG -gcp  0 0 -86517.10966018954 443.645829377676 -gcp 2310 0 14978.910498662082 30.51726815643259 -gcp 2310 1467 15088.38925912324  -64653.61709593148 -gcp 0 1467 -87149.5276955563 -64243.31846699882 -a_ullr -86517.10966018954 443.645829377676 15088.38925912324  -64653.61709593148 tokyo10m.png tokyo10m-epsg2451.png
Input file size is 2310, 1467
0...10...20...30...40...50...60...70...80...90...100 - done.
user@debian7-vmw:~/mapdata/tokyo$ gdalinfo tokyo10m-epsg2451.png
Driver: PNG/Portable Network Graphics
Files: tokyo10m-epsg2451.png
       tokyo10m-epsg2451.png.aux.xml
Size is 2310, 1467
Coordinate System is `'
Origin = (-86517.109660189540591,443.645829377676023)
Pixel Size = (43.985064467234970,-44.374412355357293)
GCP Projection =
GCP[  0]: Id=, Info=
          (0,0) -> (-86517.10966019,443.6458293777,0)
GCP[  1]: Id=, Info=
          (2310,0) -> (14978.91049866,30.51726815643,0)
GCP[  2]: Id=, Info=
          (2310,1467) -> (15088.38925912,-64653.61709593,0)
GCP[  3]: Id=, Info=
          (0,1467) -> (-87149.52769556,-64243.318467,0)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  -86517.110,     443.646)
Lower Left  (  -86517.110,  -64653.617)
Upper Right (   15088.389,     443.646)
Lower Right (   15088.389,  -64653.617)
Center      (  -35714.360,  -32104.986)
Band 1 Block=2310x1 Type=Byte, ColorInterp=Red
Band 2 Block=2310x1 Type=Byte, ColorInterp=Green
Band 3 Block=2310x1 Type=Byte, ColorInterp=Blue

この画像ファイルでも地図表示はできますが、投影法が確認できるように追加します。

user@debian7-vmw:~/mapdata/tokyo$ gdal_translate -of PNG -a_srs "EPSG:2451" tokyo10m-epsg2451.png tokyo10m-epsg2451-2.png
Input file size is 2310, 1467
0...10...20...30...40...50...60...70...80...90...100 - done.
user@debian7-vmw:~/mapdata/tokyo$ gdalinfo tokyo10m-epsg2451-2.png
Driver: PNG/Portable Network Graphics
Files: tokyo10m-epsg2451-2.png
       tokyo10m-epsg2451-2.png.aux.xml
Size is 2310, 1467
Coordinate System is:
PROJCS["JGD2000 / Japan Plane Rectangular CS IX",
    GEOGCS["JGD2000",
        DATUM["Japanese_Geodetic_Datum_2000",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6612"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4612"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",36],
    PARAMETER["central_meridian",139.8333333333333],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",NORTH],
    AXIS["Y",EAST],
    AUTHORITY["EPSG","2451"]]
Origin = (-86517.109660189540591,443.645829377676023)
Pixel Size = (43.985064467234970,-44.374412355357293)
GCP Projection =
GCP[  0]: Id=, Info=
          (0,0) -> (-86517.10966019,443.6458293777,0)
GCP[  1]: Id=, Info=
          (2310,0) -> (14978.91049866,30.51726815643,0)
GCP[  2]: Id=, Info=
          (2310,1467) -> (15088.38925912,-64653.61709593,0)
GCP[  3]: Id=, Info=
          (0,1467) -> (-87149.52769556,-64243.318467,0)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  -86517.110,     443.646) (138d52'25.30"E, 36d 0' 0.58"N)
Lower Left  (  -86517.110,  -64653.617) (138d52'50.52"E, 35d24'48.51"N)
Upper Right (   15088.389,     443.646) (140d 0' 2.53"E, 36d 0'13.97"N)
Lower Right (   15088.389,  -64653.617) (139d59'58.13"E, 35d25' 1.63"N)
Center      (  -35714.360,  -32104.986) (139d26'19.06"E, 35d42'35.91"N)
Band 1 Block=2310x1 Type=Byte, ColorInterp=Red
Band 2 Block=2310x1 Type=Byte, ColorInterp=Green
Band 3 Block=2310x1 Type=Byte, ColorInterp=Blue

nippon_bmi_akiruno_pgis.map に 次のようにレイヤ（tokyo_height）を設定します。

---
LAYER
 NAME tokyo_height
 TYPE RASTER
 STATUS ON
 DATA "../mapdata/tokyo/tokyo10m-epsg2451-2.png"
 MINSCALEDENOM 1000    # 不適当な縮尺で使用されないように。
 MAXSCALEDENOM 1000000
 METADATA
  "wms_title" "Tokyo Height 10m WMS LAYER"
  "wms_srs" "EPSG:2451"
 END
END

---

ol003-nippon_bmi_akiruno.pgis.html に次のように 東京都の高度数値モデルの画像レイヤ（tokyo_height）を layer0 に設定します。

---
   layer0 = new OpenLayers.Layer.WMS( "Tokyo Height WMS",
    "http://192.168.1.200/cgi-bin/mapserv?",
    {
     map: '/home/user/mapfile/nippon_bmi_tokyo_pgis.map',
     layers: 'tokyo_height',
     format: 'image/png'
   });

---

地図を拡大するため、コントロールパネルの "+" を4回クリックすると、公共施設の位置が表示されました。
nippon_nlni_tokyo_pgis.map の tokyo_pf レイヤと nippon_bmi_akiruno_pgis.map の akiruno_kenchiku レイヤの表示倍率を揃えるために MINSCALEDENOM と MAXSCALEDENOM の値を同じにします。（MINSCALEDENOM 500、MAXSCALEDENOM 10000 にしました。）
投稿者 時刻:
ラベル: , , , ,

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)