国土地理院データで遊ぶその2

この記事は以下の記事の続きです．

前回は，緯度経度を指定して，任意のズームレベルにてその場所の地図画像を一枚取得するところまで作成しました．

今回は，前回取得したデータと，別のデータを合わせて表示してみます．

国土地理院のサイトでは、標準地図タイルだけでなく，様々な情報が公開されています．

その中の一つに，「標高タイル」というものがあります．

標高タイル（基盤地図情報数値標高モデル）
URL：https://cyberjapandata.gsi.go.jp/xyz/dem/{z}/{x}/{y}.txt（DEM10B テキスト形式）
https://maps.gsi.go.jp/development/ichiran.html#dem

zxyの仕様は，前回取得したマップタイルのものと全く同様であるように見えます．以下のページに詳細仕様が記述されていました．

地理院地図｜標高タイルの詳細仕様

256ピクセル×256ピクセルの地図タイルと関連づけるために、標高タイルは以下のデータ仕様になっています。
1行にカンマ区切りで256個の標高値を表す数値データが格納されている。
1枚の標高タイルは256行からなる。
数値データは対応する地図タイルのピクセル座標における標高値を表す。座標系はタイル座標と同一である。つまり、東方向がX、南方向がYになっている。
標高値が存在しない画素には「e」の文字が格納されている。
標高データは小数点第二位までデータとして入っている（単位はm）。
https://maps.gsi.go.jp/development/demtile.html

つまり，マップタイル１ピクセルごとに，対応する標高情報が取得できるようです．

ということは，この標高情報を元に3Dレリーフモデルを作成し，それに地図画像をテクスチャで貼り付けてやると，簡単に立体地図が作成できるはずです．早速試してみましょう．

まず、標高が判りやすい地域を適当に選びます。判りやすい例として，富士山が含まれるタイルにしてみます．

富士山の座標は

lat : 35.3606247
long : 138.7186086

とのことです．これを前回の式で計算してみます．今回は範囲を広めに，ズームレベル10で計算します．

x : 906.5773755733335
y : 404.34916098495285

つまりこの例では、z=10, x=906, y=404 になります．実際に取得した画像は以下です．

https://cyberjapandata.gsi.go.jp/xyz/std/10/906/404.png

この範囲の標高パネルを取得するURLは以下です．

https://cyberjapandata.gsi.go.jp/xyz/dem/10/906/404.txt

取得したデータには，一行あたりカンマで区切られた256個の標高データが256行記述されています．

このCSVデータを読み込み，以下のような処理で3Dモデルを作成，それにテクスチャを貼り付けます．（詳細は省略しますがコードを読めば判ると思います）

// xにX番号，yにY番号、zoomにズームレベル，
// elavationに座標情報(256×256個)が格納されている
function drawTile(x, y, zoom, elevation){
  //地図画像ロード
  var img = "http://cyberjapandata.gsi.go.jp/xyz/std/" + zoom + "/" + x + "/" + y + ".png";
  THREE.ImageUtils.crossOrigin = "*";     //他サーバのイメージを読み込む
  var texture = THREE.ImageUtils.loadTexture(img);
  var material = new THREE.MeshLambertMaterial({map: texture});
  var gw = elevation[0].length - 1;
  var gh = elevation.length - 1; //単位高さを計算
  // 地球の赤道周りは約4万km、それを2をズームレベルで乗じたもの×(タイル一枚のピクセル数)で割って1px当たりの高さ(m)を算出
  var h = 1 / (((40000 * 1000) / (Math.pow(2, zoom) * tileW)));
  //ジオメトリを作成
  var geometry = new THREE.PlaneBufferGeometry(tileW, tileH, gw, gh);
  var vertices = geometry.attributes.position.array;
  for(var y = 0; y < (gh+1); y++) {
    for(var x = 0; x < (gw+1); x++){
      //標高データをセット
      if(elevation[y][x] != "e") {
        vertices[(y * (gw+1) + x) * 3 + 2] = parseFloat(elevation[y][x]) * h;
      } else {
         vertices[(y * (gw+1) + x) * 3 + 2] = -5 * h;
      }
    }
  }
  //法線ベクトルの自動計算
  geometry.computeFaceNormals();
  geometry.computeVertexNormals();

  //メッシュを作ってシーンに追加
  var mesh = new THREE.Mesh(geometry, material);
  scene.add(mesh);
}

その他の処理を一通り実装し，表示できるようにしたものが以下です．