USGS(土地利用データ)の読み方(FORTRAN文です)
integer lai,lao,loi,loo
parameter(ni=43200, nj=21600)
parameter(lai=153000, lao=138024)
parameter(loi=347400, loo=367236)
integer ulxmap, ulymap, xdim, ydim
integer i, j, a, b, c, d, u, v, w, z
integer*1 data(ni)
real x(ni), y(nj)
integer land(ni)
ulxmap = -647985
ulymap = 323985
xdim = 30
ydim = 30
u= (loi - ulxmap) / xdim + 1
v= (loo - ulxmap) / xdim + 1
w= (ulymap - lai) / ydim + 1
z= (ulymap - lao) / ydim + 1
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 file open
open( unit=10, form='unformatted', access='direct', recl=ni*1,
+ file='gusgs2_0ll.img', err=999 )
open( unit=20, form='unformatted', access='direct',
+ recl=((v-u) + 1) *4, file='land1.bin',err=999 )
C file read & write and dataout put
do j=z,w,-1
read( unit=10, rec=j )( data(i), i=1,ni)
do i=u,v
land(i) = data(i)
if(land(i).eq.1.or.land(i).eq.7) then
land(i)=19
else if(land(i).eq.3.or.land(i).eq.4) then
land(i)=19
else if(land(i).eq.2.or.land(i).eq.5) then
land(i)=19
else if(land(i).eq.6.or.land(i).eq.8.or.land(i).eq.9
+ .or.land(i).eq.11.or.land(i).eq.12
+ .or.land(i).eq.13.or.land(i).eq.14
+ .or.land(i).eq.15.or.land(i).eq.18) then
land(i)=19
else if(land(i).eq.10) then
land(i)=19
else if(land(i).eq.16) then
land(i)=-1
else if(land(i).eq.17) then
land(i)=19
else if(land(i).eq.19) then
land(i)=39
else if(land(i).eq.20.or.land(i).eq.21.or.
+ land(i).eq.22.or.land(i).eq.23) then
land(i)=39
else if(land(i).eq.24) then
land(i)=-1
else if(land(i).eq.99) then
land(i)=90
else if(land(i).eq.100) then
land(i)=91
else
write(*,*) 'miss'
end if
end do
do i=u,v
if(i.gt.33601.and.i.lt.33691.and.j.gt.6090.and.
+ j.lt.6156) then
land(i)=19
else
land(i)=39
end if
end do
do i=u,v
write( unit=20, rec=z-j+1)(land(i), i=u,v)
end do
do i=u,v
write(*,*) x(i)/3600, y(j)/3600, land(i)
C write(*,*) land(i)
end do
end do
close(20)
close(10)
C write(*,*) u,v,w,z
stop
999 write(*,*) 'can not open'
end