日本の斜度図

by 宮城磯治                                  平成１４年３月１８日(月)


国土地理院のデジタル標高データ(DEM 50m mesh)から、斜度データを作成しま
した。個々の斜度図ファイルは、20万分の一区画毎にとりまとめて(tarによる)、
圧縮(gzipによる)してあります。斜度図のファイルフォーマットはDEMと同様
で、標高の代わりに斜度が入力されています。その値[メートル相当] は斜度
[度] ×40です。すなわち斜度図ファイルをDEMに見たてた際、0度が標高0メー
トル、斜度90度が標高3600m に相当します。

各ピクセルの斜度算出方法は、perlサブルーチン`geisha'によりました。斜度
算出の際に用いる水平距離には、DEMのピクセルサイズは緯度によって異なる
ため、サブルーチン`mapRatio' を用いて算出した値を使いました。


# サブルーチン「geisha」; DEMの素子の傾きを計算する。
# 入力は以下の配置の標高データと、
# 素子の東西・南北サイズ[m] (geishaに依存)
# A B C
# D E F
# G H I
# 出力[m]はEの素子の傾き[deg]の40倍。
# 参考文献：http://isweb20.infoseek.co.jp/area/okinawaw/okinawa-u/ronbun4/GIS99_9.htm
sub geisha(){
  my($A,$B,$C,$D,$E,$F,$G,$H,$I,$NS,$EW) = @_;
  my($Z1,$Z2,$deg);
  $Z1= (($A + 2.0*$D + $G) - ($C + 2.0*$F + $I)) / 4.0 / 10.0;
  $Z2= (($G + 2.0*$H + $I) - ($A + 2.0*$B + $C)) / 4.0 / 10.0;
  $X1= $EW * 2.0;
  $Y1= $NS * 2.0;
  $deg= ($Z1/$X1)*($Z1/$X1) + ($Z2/$Y1)*($Z2/$Y1);
  $deg= sqrt($deg);
  $deg= atan2(1.0, $deg);  # atan2(Y,X) Returns the arctangent of 
                           # Y/X in the range -pi to pi.
  $deg= $deg * 180.0 / 3.14159265358979;
  $deg= 90.0 - $deg;
  $deg= $deg * 400.0;      #DEMは0.1m単位なので、40倍に相当。
  return($deg);
}


# サブルーチン「mapRatio」; DEMの素子の大きさをgeishaに与える。
# 入力=緯度 [deg]
# 出力=その緯度での、DEM 1pixelの距離(南北1.5秒[m],東西2.25秒[m])
# 参考文献：地球楕円体に関する計算式(理科年表、1995年、621ページ)
sub mapRatio () {
    my($Ido) = @_;
    my($L1, $L2, $L4, $Ratio);
    my($Pi, $Rad, $E2) = (3.14159265358979, 6378137, 0.006694470);

    $Ido /= 180.0;
    $Ido *= $Pi;
    
    #東西方向2.25秒の距離(m)
    $L1 = $Pi / 180000.0;
    $L1 *= $Rad * cos($Ido);
    $L1 /= sqrt(1.0-($E2 * (sin($Ido))*(sin($Ido)) ));
    $L1 /= 3600.0;
    $L1 *= 2.25;
    $L1 *= 1000.0;
    
    #南北方向1.5秒の距離(m)
    $L4 = $Pi / 180000.0;
    $L4 *= $Rad * (1.0-$E2);
    $L4 /= sqrt((1.0 - $E2 * (sin($Ido))*(sin($Ido)))*
                (1.0 - $E2 * (sin($Ido))*(sin($Ido)))*
                (1.0 - $E2 * (sin($Ido))*(sin($Ido))) );
    $L4 /= 3600.0;
    $L4 *= 1.5;
    $L4 *= 1000.0;
    
    $Ratio = $L1/$L4;
    
    return($L4, $L1, $Ratio);
    #L4が南北(ｍ/2.25秒)、L1が東西(ｍ/1.5秒)
}