[参考] DEM40 ファイル群の生成方法


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