[Info] How to generate DEM40 data files


Among program group GDMP, program gtopo creates grid data of topographic height from DEM40 data files. Here, DEM40 is a group of files of terrain elevation data in Japan, reconstructed from Digital Maps (50m mesh elevation, Japan-I, Japan-II and Japan-III) by the Geographical Survey Institute, Japan (GSI). The following process is applied in the reconstruction.

  1. Geodetic system conversion from Tokyo Datum to WGS84, including the correction for the distortion of GSI's Tokyo Datum,
  2. Mesh interpolation from [1/40 min.Lat. × 3/80 min.Long.] into [1/40 min.lat. × 1/40 min.Long],
  3. Conversion from each 'standard area mesh' [40 min.Lat. × 1 deg.Long.] file to each [1 deg.Lat. × 1 deg.Long.] area file,
  4. Conversion from the center value of mesh-element to the value at grid lines crossing,
  5. Replacing negative land elevation to 0., and sea area data to -1.,
  6. Generated data files are divided into subdirectories of UTM zone.

The source code "genDEM40m.f90" of the Fortran program realizing the process above is presented at the bottom of this page.
  The correction for the distortion of GSI's Tokyo datum is known in a form of "Enhance_Par" provided by GSI,
      http://www.gsi.go.jp/MAP/CD-ROM/sekaitaiou/2500Conv/Enhance_Par.exe
(Self-extracting archive), which could be approximated into the corrections onto the geodetic translation by 'xw84t' subprogram (entry name 'xtw84') with deviding into several zones.
  The path names of input data to this program are given as 'GSI50dem/xxxx/xxxxyy.MEM' for the area of primary mesh code xxxx, and secondary mesh code yy. Output data file for [1 deg.Lat. × 1 deg.Long.] range with SW corner at (iideg.N, kkkdeg.E) becomes the path name of "dem40/Zjj/NiiEkkk.alt", where jj is the UTM zone nunber. Before running this program, output directory dem40 and subdirectories dem40/Z51/, dem40/Z52/, dem40/Z53/, dem40/Z54/, dem40/Z55/, dem40/Z56/ must be created, if necessary. Run time of the program would be of the order of tens minutes because of a large size storage data handling.

To prepare for using the program gtopo, move the directory dem40 above (consisting of 120 files, about 45 MB each) under the absolute directory path of '/home/SHARE/data/DEM/', as defined as constant parameter DEMDIR in the 'gtopo.f90' source.


!---- Source program of generating DEM40 data files
!
  integer :: mh(200,200)
  real :: h(2800,2000), a(2401,2401)
  character(24) :: fnam1, fnam2;
  character(8) :: area;  character(6) :: cmesh, code
  logical :: exist

  do lli=20,45;  do llk=122,153
    llj = llk/6 + 31
    llc = lli*1000 + llk
    if ((lli < 24)  .and. (llc /= 20136)) cycle
    if ((llk > 145) .and. (llc /= 24153)) cycle

    write(fnam2,'(a7,i2,a2,i2,a1,i3,a4)')   &
          'dem40/Z', llj, '/N', lli, 'E', llk, '.alt'
    kc = 0
    if ((llc == 25131) .or. (llc == 32139)) kc = 1
    if ((llc >= 24122) .and. (llc <= 24125)) kc = 1
    do i=1,2800; do k=1,2000;  h(i,k) = -1.;  enddo; enddo
    key = 0
    ins = lli*12 - 1
    kns = (llk-100)*8 - 1
    do in=ins,ins+13;  do 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 /= fnam1(15:20)) call abendm('mesh-code error')
        do n=1,200
          read(1,'(a6,i3,200i5)') code, nn, (mh(k,201-n),k=1,200)
          if (code /= cmesh) call abendm('mesh-code conflict')
          if (nn /= n) call abendm('format error')
        enddo
        close(1)
        is = (in-ins) * 200
        ks = (kn-kns) * 200
        do i=1,200;  do k=1,200
          if (mh(k,i) == -9999) cycle
          if (mh(k,i) <= 0) then
            h(is+i,ks+k) = 0.
          else
            h(is+i,ks+k) = float(mh(k,i)) / 10.
          endif
        enddo;  enddo
        key = 1
      endif
    enddo;  enddo
    if (key == 0) cycle

    key = 0
    wlat0 = float(lli*60);  tlat0 = wlat0 - 5.0125
    wlon0 = float(llk*60);  tlon0 = wlon0 - 7.51875
    do kw=1,2401;  do 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 == 1) then
        if      ((tlat >= 1920.) .and. (tlat <= 1960.) .and.   &
                 (tlon >= 8340.) .and. (tlon <= 8400.)) then
          tlat = tlat - 0.0404
          tlon = tlon - 0.0023
        else if ((tlat >= 1520.) .and. (tlat <= 1560.) .and.   &
                 (tlon >= 7860.) .and. (tlon <= 7920.)) then
          tlat = tlat + 0.2016
          tlon = tlon - 0.3138
        else if ((tlat >= 1440.) .and. (tlat <= 1500.) .and.   &
                 (tlon >= 7320.) .and. (tlon <= 7560.)) then
          if (tlon >= 7500.) then
            tlat = tlat + 0.0310
            tlon = tlon - 0.0401
          else if (tlon >= 7470.) then
            tlat = tlat + 0.1527
            tlon = tlon - 0.2867
          else if (tlon >= 7443.) then
            tlat = tlat - 0.0784
            tlon = tlon - 0.1221
          else if (tlon <= 7410.) then
            tlat = tlat - 0.0782
            tlon = tlon - 0.1196
          endif
        endif
      endif
      dlat = tlat - tlat0;  ip = ifix(dlat / 0.025)
      dlon = tlon - tlon0;  kp = ifix(dlon / 0.0375)
      if ((ip <= 0) .or. (ip >= 2800) .or.   &
          (kp <= 0) .or. (kp >= 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
      flat = dlat - float(ip)*0.025
      flon = dlon - float(kp)*0.0375
      if ((h(ip,kp) < 0.) .or. (h(ip+1,kp+1) < 0.) .or.   &
          (h(ip,kp+1) < 0.) .or. (h(ip+1,kp) < 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
    enddo;  enddo
    if (key == 1) then
      write(6,*) fnam2
      area = fnam2(11:17) // ' '
      open(2,file=fnam2,status='new')
      write(2,'(a8,4x,i4,4i8)') area, 199, 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 kw=1,2401
        write(2,'((f7.1,9(1x,f7.1)))') (a(iw,kw),iw=1,2401)
      enddo
      close(2)
    endif
  enddo;  enddo
  stop
  end