gtopo30の読み方(FORTRAN文です) real lai,lao,loi,loo
parameter(ni=4800, nj=6000)
parameter(lao=38.3400, lai=39.99583333333333)
parameter(loi=100.00416666666666, loo=102.010)
real ulxmap, ulymap, xdim, ydim
integer i, j, a, b, c, d, u, v, w, z
integer*2 data(ni,nj), idata
real height(ni, nj), x(ni), y(nj)
character*2 cdata, ctmp
equivalence (idata, cdata)
ulxmap = 100.00416666666666
ulymap = 39.99583333333333
xdim = 0.00833333333333
ydim = 0.00833333333333
C read DEM data
open( unit=10, form='unformatted', access='direct', recl=ni*nj*2,
+ file='E100N40.DEM', err=999 )
read(unit=10, rec=1) ((data(i,j), i=1,ni), j=1,nj)
close(10)
C make lon. & lat.
do i=1,ni
x(i) = ulxmap + xdim*(i-1)
end do
do j=1,nj
y(j) = ulymap - ydim*(j-1)
end do
C byte swapping
do i=1,ni
do j=1,nj
idata = data(i,j)
ctmp(1:1) = cdata(1:1)
cdata(1:1) = cdata(2:2)
cdata(2:2) = ctmp(1:1)
height(i,j) = idata
end do
end do
C output
u= (loi - ulxmap) / xdim + 1
v= (loo - ulxmap) / xdim + 1
w= (ulymap - lai) / ydim + 1
z= (ulymap - lao) / ydim + 1
C write(*,*)u,v,w,z
do i=u,v
do j=w,z
write(*,*) x(i), y(j), height(i,j)
end do
end do
open( unit=20, form='unformatted', access='direct',
+ recl=((v-u) + 1)*((z-w) + 1)*4, file='chE100N40.DEM', err=999 )
write(unit=20, rec=1) ((height(i,j), i=u,v), j=w,z)
close(20)
stop
999 write(*,*) 'can not open'
end