弘前大学・理工学部 飯倉 善和


地理情報の基本的な処理

国土地理院では、日本全国をメッシュに分割して地図を発行しています。一番大きなものが1次メッシュで、縮尺20万分の1の地図に使われます。1次メッシュを縦(緯度方向)と横(経度方向)に8分割したものがで2次メッシュで、2万5千分の1の地図に使われています。

これらの地図はUTM座標系を用いて紙に投影されますが、、数値データとしては緯度と経度が基本となります。すなわち、ベクターデータは正規化座標が用いられますし、ラスターデータでは全体を縦と横に等分割したものが画素となります。

一般に衛星画像はUTM座標系に投影されますから、数値地図と衛星画像を重ね合わせるには緯度経度が基本の数値地図をUTM座標系に投影する必要があります。UTM座標は直交座標系で単位がメートルで表現されますから、距離はもちろん斜面の斜度や方位の計算も容易です。

  1. 数値地図50mメッシュ(標高)の読み込みと表示

    2次メッシュを縦横200等分した画素の中心の標高が一つのデータになります。標高は5桁の整数として0.1mまで表されています。例えば1234というデータであれば、その地点の標高は123.4mということです。海域は−9999で表されます。

    データはテキストデータですから、適当なエディタなどで読むことができます。1レコードが1009バイトの固定長で、全部で201レコードからなることが分かります。

    最初のレコードにはメッシュ番号や地図の名前などのヘッダ情報が書かれています。また各レコードの先頭の9バイトにもメッシュ番号が書かれています。残りの1000バイトが200地点の標高のデータです。

    このデータをIDLで読むには、以下を実行します。
    >dem=intarr(200,200) ; 標高を記録する配列(整数)を確保
    >openr,1, '614052.MEM' ;ファイル(青森県小泊)を開く
    >temp=intarr(200) ;1レコードのデータを読み込む配列を確保
    >readf,1,xxx,format='(i6)' ; 最初のレコードのメッシュ番号だけを読む。
    >for i=0,199 do begin & $ ; 残りの200レコードを読む。
    >readf,1,xxx,yyy,temp,format='(i6,i3,200i5)' & $ ; yyyはライン番号
    >image[*,i]=temp
    >close,1 ; ファイルを閉じる。

    配列の中身をチェックすれば、うまく読み込めたかどうかわかりますが
    >print,dem[*,0]

    画像にして表示することも簡単です。
    >window,0,xsize=200,ysize=200
    >!order=1
    >tvscl,dem

    ファイルを指定してデータの読み込む関数rmemをあらかじめ作っておけば、
    >dem=rmem('614052.MEM')
    と一行で済みます。

  2. JMCマップ(海岸線)の読み込みと表示

    JMCマップは固定長(72バイト)のテキストデータで一つのファイルに1次メッシュの行政界、道路、鉄道、河川などの情報がおさめられています。ファイル内では、2次メッシュ単位にデータが整理されており、それぞれの先頭にメッシュヘッダレコードが必ずあります。小泊の場合には以下のようになります。

    M 614052小泊                  4   10   16    8    3  183 
    
    ベクターデータは一般にノードN(端点)、ラインL(線)、エリアA(面)から構成されます。とくに行政界などの面データはDLGというフォーマットに準拠して、これらの空間的な関係を効率的に記録されています。

    行政界のデータは後で検討することとして、ここでは海岸線のラインデータだけに着目します。ラインデータレコードはLという文字で始まります。

    L  1 5    1     0    20    20 238699999    15                           
     2491 6613 2477 6594 2479 6564 2477 6548 2463 6550 2451 6558 2446 6558  
     2436 6545 2456 6507 2479 6493 2506 6496 2525 6512 2527 6537 2503 6607  
     2491 6613                                                              
    
    Lの次の最初の数字(1)はライン種類(海岸線を含む行政界)を表します。海岸線は次の数字(5)で表されます。最後の15はこのラインが15の点から構成されるのをしめしています。次の3行が各点の対応する2次メッシュ上の正規化座標x、yで、それぞれ5桁の整数で表されています。全部で30のデータがあることを確かめて下さい。

    小泊の最初の海岸線データを読み込むには、はじめにファイルを開いて小泊のデータがある場所を捜さなければなりません。以下を実行します。

    >all=string(72)
    >type=string(1)
    >flag=1
    >openr,1,'KS6140x.DAT'
    >while flag do begin & $
    >readf,1,all & $
    >reads,all,type,num,format='(a1,1x,i6)' & $
    >if(type eq 'M') and (num eq 614052) then flag=0
    >print,type,num
    
    海岸線は最初のレイヤーのラインデータですから
    >flag=1
    >while flag do begin & $
    >readf,1,all,format='(a72)' & $
    >reads,all,type,num1,num2,count,format='(a1,2x,2i2,35x,i5)' & $
    >if(type eq 'L') and (num2 eq 5) then flag=0
    >print,all
    
    次の行からのデータを30個読み込みます。
    >line1=intarr(30)
    >readf,1,line1
    >print,line1
    
    x,y成分に分けて表示します。
    >line1=reform(line1,2,15)
    >print,line1
    >plot,line1[0,*],line1[1,*]
    
    縮尺を正しくするために
    >!x.range=[0,10000]
    >!y.range=[0,10000]
    >plot,line1[0,*],line1[1,*]
    
    小さな島を表していることが分かります。ラインデータを見つけて、海岸線を読み込む作業を続ければ海岸線をすべて表示することができます。小泊の場合は海岸線のデータは5個ありました。

    以上の操作を行う関数rksを作成しました。ファイル名、2次メッシュ番号、ラインの種別(行政界1)と2(海岸線9)を指定します。

    >lline=rks('KS6140x.DAT',614052L,1,9)
    M 614052小泊                  4   10   16    8    3  183                
    L  1 5    1     0    20    20 238699999    15                           
    L  1 5    2     0    30    30 238699999     9                           
    L  1 5    3     0    40    40 238699999    15                           
    L  1 5    4     0    50    50 238699999    14                           
    L  1 5    5     0    60    11 238599999    59                           
    L  1 5    6     0    71    60 238699999   463                           
    *number of line :        6
    5
    
    読み込んだラインの数5が出力されます。llineはポインターの配列です。個別に取り出すには
    >line=*lline[5]
    とします。道路や河川も同じように処理できます。

  3. ベクターデータとラスターデータの重ね合わせ

    DEMを表示してから、

    >dem=rmem('614052.MEM')
    >window,0,xsize=200,ysize=200
    >!order=1
    >tvscl,dem
    
    海岸線を重ねます。
    >!x.range=[0,0]
    >!y.range=[0,0]
    >plot,[0,10000],[0,10000],/noerase,/nodata,xstyle=5,ystyle=5,pos=[0,0,10000,10000]
    >for i=0,5 do begin line=*lline[i] & oplot,line[0,*],line[1,*]
    
    DEMを拡大表示してから、海岸線を重ねれば見易くなります。

  4. 緯度・経度からUTM座標座標系への相互変換

    投影変換を行うためには、地球の形を定義する必要があります。地球の形は回転楕円体で近似できますが、国によってそのパラメータが多少異なります。日本はベッセルによって示された値に準拠して測量法(第11条)で、赤道半径aと扁平率fが以下のように定義されている。

      赤道半径a:6,377,397.155 m
      扁平率f:1/299.152813

    IDLでの解析では、計算に必要なパラメータを大域変数としてはじめに定義います。

    common pdata,a,b,ee,p,aee,eee
    common adata,a0,a1,a2,a3,a4,a5 ; for ll2utm (slat)
    common bdata,b0,b1,b2,b3,b4,b5 ; for utm2ll (tlat)
    a=6377397.155D
    ee=0.006674372231315D
    p=180.0D/!DPI
    aee=a*(1.0-ee)
    eee=ee/(1.0-ee)
    a0=1.005037306048555D
    a1=0.005047849240300D
    a2=0.000010563786831D
    a3=0.000000020633322D
    a4=0.000000000038853D
    a5=0.000000000000070D
    b0=1.005037306045577D
    b1=0.002511273240647D
    b2=0.000003678785849D
    b3=0.000000007380969D
    b4=0.000000000016832D
    b5=0.000000000000041D
    
    投影変換のプログラム群を作成しました。投影変換の詳しい説明は文献に譲りますが、緯度と赤道からの距離の関係を求めることが重要です。これを近似関数(slatとtlat)で行うために必要なパラメータがa0〜a5およびb0〜b5です。

    例として緯度lat(39.0)を与えて距離distを計算します。
    >print,slat(39.2345)
    4344102.3
    逆変換すれば
    >print,tlat(4344102.3)
    39.234495
    出力での丸めをなくしても
    >print,tlat(slat(39.2345))
    39.234497
    これは単精度の結果ですが、倍精度では
    >print,tlat(slat(39.2345D)), format='(f13.9)'
    39.234500000
    となり、小数点以下9桁まで正しい答えが出ます。逆をやってみるとわかりますが、この計算は1ミリメートル以上の精度が出ます。

    緯度/経度からUTM座標(x、y)への変換にはll2utmを用います。

    > print,ll2utm([39.0D,141.0D]),format='(f12.3)'
      500000.000
     4316344.510
    
    東経141度はUTMのゾーン54の中心ですから500000(m)の値になります。なお、ゾーン内での距離の精度が平均的に向上するように、中心経線上では縮尺が0.9996倍になっていることに注意して下さい。 逆変換をしてみれば
    IDL> print,utm2ll(ll2utm([39.0D,141.0D])),format='(f13.9)'
     39.000000000
    141.000000000
    
    となります。他の座標でも試してみて下さい。

    ここで等緯度線と等経度線をUTM座標系上に書いて表示してみます。

    window,1
    xs=200000
    xe=800000
    ys=4150000
    ye=4550000
    !x.range=[xs,xe]
    !y.range=[ys,ye]
    plot,[xs,xe,xe,xs,xs],[ys,ys,ye,ye,ys],/nodata
    for i=137,145 do begin for j=37,41 do ldraw,j,i
    
    プログラムldrawはll2utm.proに含まれています。

  5. 各種内挿法とDEMの投影変換

  6. 斜度と傾斜方位の計算と表示

戻る