全球の地面を眺める

標高

とりあえずいくつもある選択肢から

• ETOPO1 Ice Surface: Grid of Earth's surface depicting the top of the Antarctic and Greenland ice sheets.
• Grid/node-registered: cells are centered on lines of latitude and longitude (21601 by 10801 cells).
• ninary float

これをただ読み出すソフトは、例えばこんな風に書ける。(gfortran43 on FreeBSD)
filename='../Binary/etopo1_ice_g.flt' で読み込みファイルを指定している事に注意!

! ncols         21601
! nrows         10801
! xllcorner     -180.00833333334
! yllcorner     -90.008333333335
! cellsize      0.01666666667
! NODATA_value  -9999
! byteorder     LSBFIRST
! 

module params
implicit none
integer, parameter :: ncols=21601, nrows=10801
real*8, parameter   :: xllcorner=-180.00833333334, yllcorner=-90.008333333335
real*8, parameter   :: cellsize= 0.01666666667
real*4, parameter   :: nodata = -9999
end module

program readetopo
use params
implicit none
!
real*4              :: height(ncols)
real*8              :: posx, posy
integer             :: i, j
character(len=80)   :: filename
!
filename='../Binary/etopo1_ice_g.flt'

open(file=filename, unit=10, status='old', &
	access='direct', recl=4*ncols, form='unformatted')


do j=1, nrows
  read(10, rec=j) (height(i), i=1, ncols)
 do i=1, ncols
 call position( i, j, posx, posy )
! if( posy > 30 .AND. posy < 40 .and. posx > 130 .and.  posx < 140 ) then
 write(*,*) posx, posy, height(i)
!  end if
 end do
end do

close(10)

end

subroutine position( i, j, posx, posy )
use params
implicit none
integer :: i,j
real*8  :: posx, posy

posx = xllcorner+(i-1)*cellsize+cellsize/2
posy = yllcorner+(j-1)*cellsize+cellsize/2

end

module params
implicit none
integer, parameter :: ncols=21601, nrows=10801
real*8, parameter   :: xllcorner=-180.00833333334, yllcorner=-90.008333333335
real*8, parameter   :: cellsize= 0.01666666667
real*4, parameter   :: nodata = -9999
end module

program readetopo
use params
implicit none
!
real*4, allocatable              :: height(:, :)
real*8              :: posx, posy
integer             :: i, j
character(len=80)   :: filename
!
filename='../Binary/etopo1_ice_g.flt'

allocate(height(1:ncols, 1:nrows) )

open(file=filename, unit=10, status='old', &
	access='direct', recl=4*ncols, form='unformatted')


do j=1, nrows
  read(10, rec=j) (height(i,j), i=1, ncols)
 do i=1, ncols
 call position( i, j, posx, posy )
! if( posy > 30 .AND. posy < 40 .and. posx > 130 .and.  posx < 140 ) then
 write(*,*) posx, posy, height(i,j)
!  end if
 end do
end do

close(10)

end

subroutine position( i, j, posx, posy )
use params
implicit none
integer :: i,j
real*8  :: posx, posy

posx = xllcorner+(i-1)*cellsize+cellsize/2
posy = yllcorner+(j-1)*cellsize+cellsize/2

end

たとえばヒストグラム