国土地理院では、日本全国をメッシュに分割して地図を発行しています。一番大きなものが１次メッシュで、縮尺２０万分の１の地図に使われます。１次メッシュを縦（緯度方向）と横（経度方向）に８分割したものがで２次メッシュで、２万５千分の１の地図に使われています。
これらの地図はＵＴＭ座標系を用いて紙に投影されますが、、数値データとしては緯度と経度が基本となります。すなわち、ベクターデータは正規化座標が用いられますし、ラスターデータでは全体を縦と横に等分割したものが画素となります。
一般に衛星画像はＵＴＭ座標系に投影されますから、数値地図と衛星画像を重ね合わせるには緯度経度が基本の数値地図をＵＴＭ座標系に投影する必要があります。ＵＴＭ座標は直交座標系で単位がメートルで表現されますから、距離はもちろん斜面の斜度や方位の計算も容易です。
２次メッシュを縦横２００等分した画素の中心の標高が一つのデータになります。標高は５桁の整数として０．１ｍまで表されています。例えば１２３４というデータであれば、その地点の標高は１２３．４mということです。海域は−９９９９で表されます。
データはテキストデータですから、適当なエディタなどで読むことができます。１レコードが１００９バイトの固定長で、全部で２０１レコードからなることが分かります。
最初のレコードにはメッシュ番号や地図の名前などのヘッダ情報が書かれています。また各レコードの先頭の９バイトにもメッシュ番号が書かれています。残りの１０００バイトが２００地点の標高のデータです。
このデータをＩＤＬで読むには、以下を実行します。
配列の中身をチェックすれば、うまく読み込めたかどうかわかりますが
画像にして表示することも簡単です。
ファイルを指定してデータの読み込む関数rmemをあらかじめ作っておけば、
JMCマップは固定長（７２バイト）のテキストデータで一つのファイルに１次メッシュの行政界、道路、鉄道、河川などの情報がおさめられています。ファイル内では、２次メッシュ単位にデータが整理されており、それぞれの先頭にメッシュヘッダレコードが必ずあります。小泊の場合には以下のようになります。
行政界のデータは後で検討することとして、ここでは海岸線のラインデータだけに着目します。ラインデータレコードはLという文字で始まります。
小泊の最初の海岸線データを読み込むには、はじめにファイルを開いて小泊のデータがある場所を捜さなければなりません。以下を実行します。
以上の操作を行う関数rksを作成しました。ファイル名、２次メッシュ番号、ラインの種別（行政界１）と２（海岸線９）を指定します。
ＤＥＭを表示してから、
投影変換を行うためには、地球の形を定義する必要があります。地球の形は回転楕円体で近似できますが、国によってそのパラメータが多少異なります。日本はベッセルによって示された値に準拠して測量法（第１１条）で、赤道半径ａと扁平率ｆが以下のように定義されている。
赤道半径ａ：6,377,397.155 m
ＩＤＬでの解析では、計算に必要なパラメータを大域変数としてはじめに定義います。
例として緯度lat（39.0)を与えて距離distを計算します。
緯度／経度からＵＴＭ座標（x、y）への変換にはll2utmを用います。
ここで等緯度線と等経度線をＵＴＭ座標系上に書いて表示してみます。
>dem=intarr(200,200) ; 標高を記録する配列（整数）を確保
>openr,1, '614052.MEM' ;ファイル（青森県小泊）を開く
>temp=intarr(200) ;１レコードのデータを読み込む配列を確保
>readf,1,xxx,format='(i6)' ; 最初のレコードのメッシュ番号だけを読む。
>for i=0,199 do begin & $ ; 残りの２００レコードを読む。
>
>
>close,1 ; ファイルを閉じる。
>print,dem[*,0]
>window,0,xsize=200,ysize=200
>!order=1
>tvscl,dem
>dem=rmem('614052.MEM')
と一行で済みます。
M 614052小泊 4 10 16 8 3 183
ベクターデータは一般にノードＮ（端点）、ラインＬ（線）、エリアＡ（面）から構成されます。とくに行政界などの面データはＤＬＧというフォーマットに準拠して、これらの空間的な関係を効率的に記録されています。
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の次の最初の数字（１）はライン種類（海岸線を含む行政界）を表します。海岸線は次の数字（５）で表されます。最後の１５はこのラインが１５の点から構成されるのをしめしています。次の３行が各点の対応する２次メッシュ上の正規化座標x、yで、それぞれ５桁の整数で表されています。全部で３０のデータがあることを確かめて下さい。
>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
次の行からのデータを３０個読み込みます。
>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,*]
小さな島を表していることが分かります。ラインデータを見つけて、海岸線を読み込む作業を続ければ海岸線をすべて表示することができます。小泊の場合は海岸線のデータは５個ありました。
>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
読み込んだラインの数５が出力されます。llineはポインターの配列です。個別に取り出すには
>line=*lline[5]
とします。道路や河川も同じように処理できます。
>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,*]
ＤＥＭを拡大表示してから、海岸線を重ねれば見易くなります。
扁平率ｆ：1/299.152813
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です。
>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
となり、小数点以下９桁まで正しい答えが出ます。逆をやってみるとわかりますが、この計算は１ミリメートル以上の精度が出ます。
> print,ll2utm([39.0D,141.0D]),format='(f12.3)'
500000.000
4316344.510
東経１４１度はＵＴＭのゾーン５４の中心ですから５０００００（m）の値になります。なお、ゾーン内での距離の精度が平均的に向上するように、中心経線上では縮尺が０．９９９６倍になっていることに注意して下さい。
逆変換をしてみれば
IDL> print,utm2ll(ll2utm([39.0D,141.0D])),format='(f13.9)'
39.000000000
141.000000000
となります。他の座標でも試してみて下さい。
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に含まれています。