Webで地球地図データ(標高第1.1版)を使う。－(その１:BILをBMP化する)

# Webで地球地図データ(標高第1.1版)を使う。－(その１:BILをBMP化する)

国土地理院が公開している様々なデータのうち、地球地図日本のデータから、標高データ(BILデータ)を可視化してWebで利用するシステムを作る過程を説明していきます。
一から作ると色々難しいので、極力簡単に有るものを使う方向です。
(地図射影法とかワールドモデルとかジオイドとかそう言う面倒な話は省略したい。)

そもそもネットに繋げれる環境なら、国土地理院の提供するサービス使えばもっと簡単なんですが、最終的にオフラインで使うため色々手動でやってます。

**(2016/5/7追記：その1投稿後、QGISでGeoTiff化が出来るのが判ったため、その2で方針変更しています。結果的に後半に示したPowerShellを使ったBMP化は無駄になりましたので、ご注意願います。)**


|ソフト等|バージョン|備考|
|:------|:--------|:---|
|地球地図データ(標高第1.1版)|gm-jpn-el_1_1.zip|[国土地理院－地球地図日本のデータ](http://www.gsi.go.jp/kankyochiri/gm_jpn.html)よりダウンロード。|
|Windows PowerShell|4.0||

## 地図データについて

地図データには大きく２種類の形式があります。

1. 海岸線や領土の境などを示すベクタデータ。

1. 一定の領域の高度や植生などを集めたラスタデータ。

ベクタデータはそれを繋げてポリゴンなどの形にしたり、内外判定したりする必要があって一から扱うのは難しいですが、ラスタデータは単純に並べただけですので処理自体は簡単です。
（今なら、ベクタデータは地図データベースで処理するのがお勧めですね。）
[国土地理院－地球地図日本のデータ](http://www.gsi.go.jp/kankyochiri/gm_jpn.html)の地球地図日本ダウンロードにある表のうち「ラスタ」とあるのがラスタデータ、「ベクタ」とあるのがベクタデータです。
それぞれ内容やデータの格納形式が異なり、データの新旧もありますので、適切なものをダウンロードします。

## 今回使うデータ

今回は[標高データの第1.1版ラスタのBILデータ](http://www1.gsi.go.jp/geowww/globalmap-gsi/download/data/gm-japan/gm-jpn-el_1_1.zip)を使います。
早速解凍すると、以下のような構成になっています。

> gm-jpn-el_1_1
> └ globalm
> 　└ area
> 　　└ raster
> 　　　└ elwhaf.bil
> 　　　└ elwhaf.hdr

今回使うのはrasterフォルダ内の、\*.bilと\*.hdrのファイルだけです。
これらファイルに収録されている範囲は次の画像で色のついている範囲になります。
（表示には、ISCGMのGM Viewer ver.2を使っています。登録が必要ですが[ここのviewer](https://www.iscgm.org/gmd/)から入手できます。）
なお、あくまで日本だけ収録されているデータですので、国外は収録されていません。
（位置的に左下に台湾や中国の一部が含まれるはずですが、表示の通りデータ上は海になっています。海外も必要なら、[ここのGlobal Elevation](https://www.iscgm.org/gmd/)から世界中のデータが入手できますが、ファイルの精度が違うので注意してください。）
![GMV2.png](https://qiita-image-store.s3.amazonaws.com/0/111541/e8c19cd2-d07f-43ce-405f-3b97c8407307.png)

各ファイル形式の正式な仕様はISCGMで[公開](http://www.iscgm.org/gm/spec.html)されています。(英語ですが)
以下は必要なところだけ、かいつまんだ説明になります。

## HDRファイルの内容

HDRファイルはテキストなので、エディタ等で開けば見えますし、意味も大体判りやすいです。

```elwhaf.hdr
BYTEORDER M
LAYOUT BIL
NROWS 600
NCOLS 600
NBANDS 1
NBITS 16
BANDROWBYTES 1200
TOTALROWBYTES 1200
BANDGAPBYTES 0 
NODATA -9999
ULXMAP 120.00416666666
ULYMAP 24.995833333334
XDIM 0.00833333333333
YDIM 0.00833333333333
```

* BYTEORDER：高度データは、後にNBITSで示すだけ(16BIT=2BYTE)で格納されていますが、その2BYTEの順番を示しています。M=Big endianで格納されています。他の地図データには、Iとなっているものがあり、こちらはLittle edianで格納されています。
（M=Motorola、I=Intelの頭文字で、それぞれのCPUのメモリ格納方式に由来）
* NROWS、NCOLS：BILファイルの中に格納されている高度が、何行何列分あるか示しています。
* NODATA：高度データの無い場所を示す値です。基本的には-9999＝海を示すことになります。なお、内陸の水域は区別されません。（例：川とか湖とかは-9999ではなく、水面相当の高度が設定されています。区別するには、水域のベクタデータを使って判断する必要がありますが、今回は省略です。）
* ULXMAP、ULYMAP：BILファイルに含まれている高度がどこのものなのか、領域の北西端（上左端）の緯度、経度で示しています。
* XDIM、YDIM：1つの高度が、どれだけの領域を示すか、緯度角度、経度角度で示しています。

なお地球地図データ\(標高第1.1版\)では、ULXMAP、ULYMAP以外は値が変わらないため、以降の処理ではHDRファイルを参照せずに、処理を簡略化しています。
もし他のBILデータを使う場合は適切に変更してください。

## BILファイルの内容

BILファイルをバイナリエディタで開いてみました。

![sea.png](https://qiita-image-store.s3.amazonaws.com/0/111541/a4a749c6-f2e9-d143-3213-5f89505c6236.png)
 
日本は海が多いので、適当に開くと「D8 F1」が並んでいるところが多いです。
この「D8 F1」を10進数に直すと「-9999」となり、海であることが判ります。
また、このファイルの中に格納されている高度は、領域の南西から順に、以下のような順で格納されています。

![BIL.png](https://qiita-image-store.s3.amazonaws.com/0/111541/adc543ca-19e4-155e-f135-32e287320574.png)

先のHDRファイルで、ULXMAPとULYMAPで示されている座標は、図の赤点の位置になります。
また、XDIMとYDIMの関係も図の通りとなります。

なお、ここで唯一気を付けなければならないのは、XDIMやYDIMはあくまで緯度、経度の角度距離であり、単純にメートル換算できないことです。
やや古い地図データの説明で「250mメッシュ」などと言われるデータがありますが、あれはあくまで１マス（ここでいうXDIMやYDIM）が大体250mぐらいになっていると言うだけで、正確に250ｍではありません。
ごく狭い範囲をイメージ的に表示するだけなら気にすることはありませんが、ある程度の範囲でちゃんとした処理をする際には避けて通れない違いなので、頭の片隅に入れておいてください。

## BILをBMP化する

最終的にはGeoTiffとかの、扱いやすい形式にしたいのですが、今回はまずBILファイルを丸ごと単純なBMPにするところまでまとめます。
そんな難しいことはしないのでツールは何でも良いんですが、手ごろなPowerShellでやってみました。
色々手抜きですが、BILファイル名を指定するとBMPファイルを作ってくれます。

```bil2bmp.ps1
Add-Type -AssemblyName System.Drawing

$bil_fn = (Get-Item -Path $args[0]).FullName;
$bmp_fn = Join-path -Path ([System.IO.Path]::GetDirectoryName($bil_fn)) `
  -ChildPath ([System.IO.Path]::GetFileNameWithoutExtension($bil_fn) + ".bmp");

# 定数
$NROWS = 600;
$NCOLS = 600;
$TOTALROWBYTES = 1200;
$NODATA = -9999;

# 高度帯毎のColor準備
$c0 = [System.Drawing.Color]::FromArgb(0, 0, 0, 0);      # sea(=NODATA)
$c1 = [System.Drawing.Color]::FromArgb(255, 0, 255, 0);  # under 1000m
$c2 = [System.Drawing.Color]::FromArgb(255, 0, 0, 255);  # 1001-2000m
$c99 = [System.Drawing.Color]::FromArgb(255, 255, 0, 0);  # other

# BIL読み込み
write-host "# $bil_fn を読み込み中";
$bil = Get-Content -Path $bil_fn -Encoding Byte -ReadCount $TOTALROWBYTES

# BMP作成
$bmp = New-Object System.Drawing.Bitmap($NCOLS, $NROWS);
for($row=0 ; $row -lt $NROWS; $row++ ) {

  # 進捗表示
  $p = [int]($row / $NROWS * 100);
  Write-Progress "BMP作成" ([String]$p + "%") -percentComplete $p

  for($col=0 ; $col -lt $NCOLS; $col++ ) {

    # 高度算出
    $hex = $bil[$row][$col*2]*256 + $bil[$row][$col*2+1];
    if($hex -band 0x8000) {
      $alt = -((($hex -band 0x7fff) -bxor 0x7fff) + 1)
    } else {
      $alt = $hex;
    }

    # BMP色付け
    if($alt -eq $NODATA) {
      $bmp.SetPixel($col ,$row ,$c0);
    } elseif($alt -le 1000) {
      $bmp.SetPixel($col ,$row ,$c1);
    } elseif($alt -le 2000) {
      $bmp.SetPixel($col ,$row ,$c2);
    } else {
      $bmp.SetPixel($col ,$row ,$c99);
    }
  }
}

# BMP書き込み
write-host "# $bmp_fn を書き込み中";
$bmp.Save($bmp_fn);
write-host "# BMP書き込み終了";
```

各定数は、HDRファイルを読み込むのが理想なのですが、手を抜いてます。
高度帯毎の色はもっと細かく増やせば綺麗な画像になります。
高速化は考えてないので遅いです。プログレスバーが出るので、ちょっと待ちましょう。

## 実行結果

> PS C:\work> .\bil2bmp.ps1 .\elxjaf.bil
> \# C:\work\elxjaf.bil を読み込み中
> \# C:\work\elxjaf.bmp を書き込み中
> \# BMP書き込み終了

elxjaf.bmpをそのままpng化したのが次の画像です。

![elxjaf.png](https://qiita-image-store.s3.amazonaws.com/0/111541/e7774834-6487-6d68-b26f-13085e4beacc.png)

高度帯が大ざっぱで、色も適当ですが、処理できていることが判ると思います。
（上で言っていたように内陸の水域は区別されていないので、琵琶湖とかも陸地扱いになっています。）
今後画像を使う方法を構築するのに合わせて、色使いや画像フォーマットは変更していきたいと思います。