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