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.
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