ヤマレコBLOG: GPSデータからグラフ作り：４）経度・緯度で指定した各地点の標高データを求める(SRTMから標高取得)

« GPSデータからグラフ作り：３）経度・緯度で指定した2地点間の距離を求める | メイン | now loadingの画像を生成してくれるAjaxload »

2007年03月11日

GPSデータからグラフ作り：４）経度・緯度で指定した各地点の標高データを求める(SRTMから標高取得)

前回まででGPXファイルの標高データをもとにグラフを描くことができました。ただ、標高データがないGPXファイルの場合はグラフを描くことができないので、別の標高データを取ってくる必要があります。

国土地理院の提供する標高情報で言うと「数値地図５０ｍメッシュ（標高）」という標高データがあります。が、これをサービスとして提供するには結構な手続きが必要になります。
具体的にはCD-ROM3枚を購入した上で、測量法第２９条又は第３０条の規定に従って「測量成果の複製・使用承認手続」が必要になります。申請の際には開発目的とか開発計画を届け出て承認が出るまで数週間待つようです。

そこでもっと簡単に使える標高データがないか探してみました。
最終的に使えそうなデータとして
・国土地理院の地球地図日本
というデータと、
・NASAの提供するスペースシャトル地形データ（SRTM）
の2つぐらいに絞れました。

どちらのデータも登山ルートの標高グラフを作ってみたのですが、地球地図の方の標高メッシュデータは精度が低すぎてグラフが本来の形とはかけ離れたものになってしまいました。

SRTMを使った標高取得プログラムを作って、グラフ作成をしました。
標高データの入手からプログラム作成までの流れをまとめると、、、

１）標高データを入手
SRTMの標高データはFTPサーバ（ftp://e0srp01u.ecs.nasa.gov/srtm/）から入手できます。

２）データフォーマットを理解
地形データは緯度経度ごとに分けられていて、例えばN37E138.hgtファイルには北緯37度、東経138度～北緯38度、東経１39度までのデータが入っています。日本のデータはEurasiaの中にあります。
データフォーマットなどの詳細は
ftp://e0srp01u.ecs.nasa.gov/srtm/version2/Documentation/SRTM_Topo.pdf
に書かれています。

仕様書とデータを見て、分かったところだけざっくり説明すると、（間違っているかもしれませんが・・・）

・各ファイルには縦に1201個、横に1201個のデータが入っていて、境界部分のデータは2つのファイルに同じデータが重なって入っている。
　⇒ 一番北西の標高データから東に向かって1201個の標高データが並んだ後、一段南に下がって次の1201個のデータが西から東に向かってならんでいくイメージです。
srtmの説明１

・データにヘッダはなく、標高データだけがバイナリ形式で16バイトのBigEndianで並んでいる。
　⇒ 例えば0x0115なら、277mという感じ。

サンプルコードも貼り付けておきます。
LatLng2Height:::緯度経度から標高算出APIのような機能を自分のサーバ内で作りたい方はご利用ください。

<html><head><title>SRTM標高計算サンプル</title></head><body>
<?php

$lat=36.3423;   // 経度：度で指定
$lon=138.1245;  // 緯度：度で指定

$b1=intval($lat);
$l1=intval($lon);
$fn = "./srtm/N".$b1."E".$l1.".hgt";

if(is_file($fn)) {
  $b=$lat;
  $l=$lon;
  $ofs = (floor((1-($b-$b1))*1200)*1201 + floor(($l-$l1)*1200)) * 2;

$fhandle = fopen($fn, "rb");
  #左上の標高
  fseek($fhandle, $ofs);
  $h1 = unpack("n", fread($fhandle, 2));
  #右上の標高
  $h2 = unpack("n", fread($fhandle, 2));
  #左下の標高
  fseek($fhandle, ($ofs+2402));
  $h3 = unpack("n", fread($fhandle, 2));
  #右下の標高
  $h4 = unpack("n", fread($fhandle, 2));

$tmp=1-($b-$b1)*1200;
  $d=$tmp-floor($tmp); #下との距離
  $c=1-$d; #上との距離
  $tmp=($l-$l1)*1200; 
  $a=$tmp-floor($tmp); #左との距離
  $b=1-$a; #右との距離

$ele1 = $a/($a+$b)*($h2[1]-$h1[1])+$h1[1];
  $ele2 = $a/($a+$b)*($h4[1]-$h3[1])+$h3[1];
  $ele = round($c/($c+$d)*($ele2-$ele1)+$ele1);

print "【座標】 \n";
 print "経度：".$lat." - 緯度：".$lon." \n \n";
 print "【四隅の標高】 \n左上：".$h1[1]."m - 右上：".$h2[1]."m \n";
 print "左下：".$h3[1]."m - 右下：".$h4[1]."m \n \n";
 print "【距離】 \n左".$a." \n";
 print "右".$b." \n";
 print "上".$c." \n";
 print "下".$d." \n \n";
 print "【結果】：".$ele."m";
} else {
 print "file open error:".$fn;
}

?>
</body></html>

サンプルコードを実行してみる。
サンプルコードをダウンロード

あとはSRTMのデータの欠落部分を修正し、メッシュデータと座標位置との変換をしつつグラフを作ることができました。

↑完成したグラフ

投稿者 matoyan : 2007年03月11日 21:13

トラックバック

このエントリーのトラックバックURL:
http://www.yamareco.com/weblog/mt-tb.cgi/10

ヤマレコBLOG

2007年03月11日

GPSデータからグラフ作り：４）経度・緯度で指定した各地点の標高データを求める(SRTMから標高取得)

トラックバック

コメント

コメントしてください

カテゴリー

最近のエントリー

アーカイブ