GDMP プログラム群中の gtopo プログラムは,DEM40 データから
地形高度グリッドデータを作成するが,その DEM40 データとは,国土地理院発行の
数値地図50mメッシュ (標高) の CD-ROM 3枚(日本-T,日本-U,日本-V)
[訂正情報を含む] を用いて,日本の地表高データを再構成したものである.
その再構成の内容は,
(1) 旧東京測地系から世界測地系へ,旧東京測地系の歪み補正を含めて 変換,
(2) Mesh size を,緯度 1/40 分,経度 3/80 分 から 緯経度とも 1/40分 に,
(3) 緯度範囲・経度範囲とも 1度 の区画ごとのデータファイルに,
(4) 細分Mesh中心位置の値 から 細分Meshの格子点での値へ,
(5) 負の陸域標高値を 0. に置換,海域を示すデータを -1. に置換,
である.
この処理を行う,FORTRANプログラムのソースコードを 下記に示す.
旧東京測地系の歪み補正については,
国土地理院が提供している「拡張パラメータファイル」
http://www.gsi.go.jp/MAP/CD-ROM/sekaitaiou/2500Conv/Enhance_Par.exe
の情報を参考に,必要な精度を考慮して,数多くない地域に区分して 測地系変換
サブプログラム 'xtw84' の結果に必要な補正を加えるようにしてある.
このプログラムの入力データは,国土地理院による元データにならって,
1次メッシュコードが 'xxxx' で 2次メッシュコードが 'yy'
で表される範囲のデータが ファイルのpath名
'GSI50dem/xxxx/xxxxyy.MEM' で表現されるものとし,
出力データは,緯度 'ii' 度・経度 'kkk' 度 を南西角とする
緯経度とも 1度 の範囲のデータが,path名で 'wm40/iikkk.alt'
のファイルとなる.(すなわち,GSJ50dem サブディレクトリからデータ入力し,
wm40 サブディレクトリへデータ出力する.)
gtopo プログラムでの利用のためには, 上記の wm40 サブディレクトリへ出力される 120個のファイル (1ファイル約45MBの固定サイズ)を,絶対path が '/home/SHARE/data/DEM/wm40/' なるディレクトリ内に配置すればよい.
c---- Source program of generating DEM40 data files c dimension mh(200,200), h(2800,2000), a(2401,2401) character fnam1*24, fnam2*14, cmesh*6, code*6, area*8 logical exist c do 100 lli=20,45 do 100 llk=122,153 llc = lli*1000 + llk write(fnam2,'(a5,i5,a4)') 'wm40/', llc, '.alt' if (lli.lt.24 .and. llc.ne.20136) goto 100 if (llk.gt.145 .and. llc.ne.24153) goto 100 c kc = 0 if (llc.eq.25131 .or. llc.eq.32139) kc = 1 if (llc.ge.24122 .and. llc.le.24125) kc = 1 do 200 i=1,2800 do 200 k=1,2000 h(i,k) = -1. 200 continue key = 0 ins = lli*12 - 1 kns = (llk-100)*8 - 1 do 300 in=ins,ins+13 do 300 kn=kns,kns+9 ima = in / 8 imb = in - ima*8 kma = kn / 8 kmb = kn - kma*8 write(fnam1,'(a9,2i2,a1,2i2,2i1,a4)') * 'GSI50dem/', ima,kma, '/', ima,kma,imb,kmb, '.MEM' inquire(file=fnam1,exist=exist) if (exist) then open(1,file=fnam1,status='old') read(1,'(a6)') cmesh if (cmesh.ne.fnam1(15:20)) call abendm('mesh-code error') do 400 n=1,200 read(1,'(a6,i3,200i5)') code, nn, (mh(k,201-n),k=1,200) if (code.ne.cmesh) call abendm('mesh-code conflict') if (nn.ne.n) call abendm('format error') 400 continue close(1) is = (in-ins) * 200 ks = (kn-kns) * 200 do 500 i=1,200 do 500 k=1,200 if (mh(k,i).eq.-9999) goto 500 if (mh(k,i).le.0) then h(is+i,ks+k) = 0. else h(is+i,ks+k) = float(mh(k,i)) / 10. endif 500 continue key = 1 endif 300 continue if (key.eq.0) goto 100 c key = 0 wlat0 = float(lli*60) wlon0 = float(llk*60) tlat0 = wlat0 - 5.0125 tlon0 = wlon0 - 7.51875 do 600 kw=1,2401 do 600 iw=1,2401 wlat = wlat0 + float(iw-1)/40. wlon = wlon0 + float(kw-1)/40. call xw84t(wlat, wlon, 0., tlat, tlon, talt) if (kc.eq.1) then if (tlat.ge.1920. .and. tlat.le.1960. .and. * tlon.ge.8340. .and. tlon.le.8400.) then tlat = tlat - 0.0404 tlon = tlon - 0.0023 else if (tlat.ge.1520. .and. tlat.le.1560. .and. * tlon.ge.7860. .and. tlon.le.7920.) then tlat = tlat + 0.2016 tlon = tlon - 0.3138 else if (tlat.ge.1440. .and. tlat.le.1500. .and. * tlon.ge.7320. .and. tlon.le.7560.) then if (tlon.ge.7500.) then tlat = tlat + 0.0310 tlon = tlon - 0.0401 else if (tlon.ge.7470.) then tlat = tlat + 0.1527 tlon = tlon - 0.2867 else if (tlon.ge.7443.) then tlat = tlat - 0.0784 tlon = tlon - 0.1221 else if (tlon.le.7410.) then tlat = tlat - 0.0782 tlon = tlon - 0.1196 endif endif endif dlat = tlat - tlat0 ip = ifix(dlat / 0.025) flat = dlat - float(ip)*0.025 dlon = tlon - tlon0 kp = ifix(dlon / 0.0375) if (ip.le.0.or.ip.ge.2800.or.kp.le.0.or.kp.ge.2000) then write(6,*) llc write(6,*) ' ip=',ip, ' kp=',kp, ' iw=',iw, ' kw=',kw write(6,*) ' tlat=',tlat, ' tlon=',tlon write(6,*) ' wlat=',wlat, ' wlon=',wlon write(6,*) ' lat0=',tlat0, ' lon0=',tlon0 call abendm('out of range') endif c flon = dlon - float(kp)*0.0375 if (h(ip,kp).lt.0..or.h(ip+1,kp+1).lt.0..or. * h(ip,kp+1).lt.0..or.h(ip+1,kp).lt.0.) then a(iw,kw) = -1. else h0 = h(ip,kp) + (h(ip+1,kp)-h(ip,kp))*flat h1 = h(ip,kp+1) + (h(ip+1,kp+1)-h(ip,kp+1))*flat a(iw,kw) = h0 + (h1-h0)*flon key = 1 endif 600 continue if (key.eq.0) goto 100 c write(6,*) fnam2 area = 'w' // fnam2(6:10) // ' ' open(2,file=fnam2,status='new') write(2,'(a8,i4,4x,4i8)') area, 399, lli*60, llk*60, 0, 0 write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') * 0,0, 25,25, 2401,2401, 9999.9, -1. do 700 kw=1,2401 write(2,'((f7.1,9(1x,f7.1)))') (a(iw,kw),iw=1,2401) 700 continue close(2) c 100 continue stop end