SRTM3データから標高を読む

NASAが公開しているSRTM3標高データファイルから，perlで値を読み込むサンプルです．

公式サイト

Shuttle Radar Topography Mission - SRTM

日本のデータファイルは，

http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/

のEurasia配下にあります．たくさんありますが，ファイル名の規則はこうなっています．

↑

ファイル名の規則

N35E139.hgt の場合

左下(南西)の隅が，北緯35度東経139度
北緯35～36度，東経139～140度の範囲を含む

となります．日本全域をカバーしようとすると，全部で270ファイルぐらいになると思います．

↑

ファイルの内容

この範囲を3秒ごとに分割した1201×1201のセルの，ジオイドモデルWGS84/EGM96での標高値が最初から最後まで並んでいるというのが，ファイルの内容です．

各標高値は，ビッグエンディアンの2バイト整数で格納されており，欠測値は-1です．

ファイルサイズはぴったり2884802(=1201×1201×2)バイトです．

↑

プログラム例

#!/usr/bin/perl

use strict;

my $COLS = 1201;
my $RECSIZE = $COLS * 2;

# 対象位置の緯度経度を，起動引数から取得する
my($y, $x) = (shift, shift);

my $fx = int($x);
my $fy = int($y);

my $filename = "N" .$fy. "E" .$fx. ".hgt";

my $base_x = $fx;
my $base_y = $fy + 1;

my $x_pos = int( ($x - $base_x) * 1200 );
my $y_pos = int( ($base_y - $y) * 1200 );

if (open(HGTFILE, $filename)) {
	binmode(HGTFILE);

	seek(HGTFILE, $RECSIZE * $y_pos + $x_pos * 2, 0);

	my $buf;
	read(HGTFILE, $buf, 2);
	close(HGTFILE);

	# ビッグエンディアンのunsigned shortで取得する
	my $h = unpack("n*", $buf);

	# unsignedからsignedに変換する
	$h = $h - 65536 if ($h > 32767);

	print "height=$h\n";
}

↑

動いているサンプル

緯度経度取得ツールで，実際に位置から標高値を取得するのを試すことができます．