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