ขั้นตอน-ปรับโครงสร้างGTOPO30 1ระวางเพื่อประมวลผลและแสดงผลใน GRADS

การจัดแต่งโครงสร้างโมเดลความสูง GTOPO30 ใน Fortran เพื่อการประมวลผลและเปิดใน GRADS

จุดประสงค์ของผมคือต้องการอ่านไฟล์ความสูง 1 ระวางใน ฟอร์แทรน เพื่อทำการประมวลผลในขั้นต่อไป แล้วต้องการเปิดใน GRADS โดยไม่ต้องการที่จะระบุ YREV ตรง Option นั่นคือ ผมต้องการจะแก้โครงสร้างการจัดเก็บข้อมูลความสูงเพื่อให้ง่ายต่อการประมวลผลในฟอรแทน

1.)   ดาวโหลดไฟล์โมเดลความสูงDEM Gtopo30 จาก https://lta.cr.usgs.gov/get_data ทำการswap dataของแต่ละTile ก่อน เช่น dd if=gt30e100n40.dem of= DEM_tile_E100N40.dat conv=swab

2.)   ใช้ฟอร์แทรนช่วยในการจัดโครงสร้างข้อมูล

program read_gtopo1tile

!-------------------------------------------------

! 2014.07.30

! read gtopo 1 tile after swap to little endian

! write to binary and plot in grads

!-------------------------------------------------

implicit none

character*150 :: infile_path

integer,parameter :: nrow=6000,ncol=4800                    ! dem dimension

integer,parameter :: row=6000,col=4800

integer,parameter :: sizeofreal=4                         !-record size

integer,parameter :: sizeofint=2                        !-record size

integer*2,dimension(ncol,nrow) :: deml,demr                    !-dem left and right

real*4,dimension(col,row) :: demn                      !-merged dem

integer :: i,j,m,n

infile_path='/home/mnattapon/RAID3storage4/research_icp/analysis_phb/beamblock_ans/dem_preparation/dem_raw_gtopo30/'

open(unit=10,file=trim(infile_path)//'E100N40.DEM2',form='unformatted',access='direct',recl=ncol*nrow*sizeofint)

print*,trim(infile_path)//'E060N40.DEM2'

read(10,rec=1)demr

close(10)

print*,'reading...'

do i=1,ncol

  do j=1,nrow

     m=ncol-i+1  !-col inverse

     n=nrow-j+1  !-row inverse

     demn(i,n)=demr(i,j)*1.0

  enddo

enddo

open(unit=20,file='DEM_tile_E100N40.dat',form='unformatted',access='direct',recl=row*col*sizeofreal)

write(20,rec=1)demn

close(20)

endprogram read_gtopo1tile

3.)   สร้าง control file ของข้อมูลเพื่อเปิดในโปรแกรม GRADS ชื่อ ctl_read_dem1tile.ctl

dset DEM_tile_E100N40.dat

title GTOPO30 tile E100N40 0.008333333333333 deg

undef -9999.0

**options little_endian

XDEF 4800  LINEAR   100     0.00833333333333

YDEF 6000  LINEAR  -10    0.00833333333333

zdef 1 levels 0 1

tdef 1 linear 22Z24aug2010 1hr

vars 1

dem     0 99 elevation [m]

endvars

4.)   เข้าโปรแกรม GRADS และเปิดข้อมูลภาพออกมาดู

open ctl_read_dem1tile.ctl

set gxout shaded

d dem

cbarn

ได้ภาพดังข้างล่างนี้

ขั้นตอน-การต่อระวางโมเดลความสูง GTOPO30 2 ระวาง

การต่อGtopo30 dem 2 tiles ในFortran และ เปิดเป็นแผนที่ดูใน GRADS
เป็นวิธีการต่อโมเดลความสูงของผมเอง โดยผมต้องการต่อเพื่อเอาไปประมวลผลความสูงกับลมพื้นผิวในขั้นต่อไปเพื่อหาฝนที่เกี่ยวเนื่องกับภูเขา โปรแกรมอาจจะไม่สวยเพราะเขียนแบบลูกทุ่ง แต่ก็ใช้งานได้ ลองเอาไปทำกันดูครับ

1.) ดาวโหลดไฟล์โมเดลความสูงDEM Gtopo30 จาก https://lta.cr.usgs.gov/get_data ทำการswap dataของแต่ละTile ก่อน เช่น   dd if=gt30e100n40.dem of= DEM_tile_E100N40.dat conv=swab
2.) ต่อ DEM ใน Fortran ดังนี้
program mosaic_dem_2tiles
implicit none
character*150 :: infile_path
integer,parameter :: nrow=6000,ncol=4800                  ! dem dimension
integer,parameter :: row=6000,col=9600                    ! dem merged dimension
integer,parameter :: sizeofreal=4                         !-record size
integer,parameter :: sizeofint=2                          !-record size
integer*2,dimension(ncol,nrow) :: deml,demr               !-dem left and right
real*4,dimension(col,row) :: demn                         !-merged dem
integer :: i,j,m,n

infile_path='/home/mnattapon/RAID3storage4/research_icp/analysis_phb/beamblock_ans/dem_preparation/dem_raw_gtopo30/'
open(unit=10,file=trim(infile_path)//'E060N40.DEM2',form='unformatted',access='direct',recl=ncol*nrow*sizeofint)
open(unit=20,file=trim(infile_path)//'E100N40.DEM2',form='unformatted',access='direct',recl=ncol*nrow*sizeofint)
read(10,rec=1)deml
read(20,rec=1)demr
close(10);close(20)
print*,'merging...'
do i=1,ncol
  do j=1,nrow
         m=ncol+i
         n=nrow-j+1  !-row inverse dem left
         demn(i,n)=deml(i,j)*1.0
         demn(m,n)=demr(i,j)*1.0
  enddo
enddo
open(unit=30,file='DEM_2tiles_merged.dat',form='unformatted',access='direct',recl=col*row*sizeofreal)
write(30,rec=1)demn
close(30)
endprogram mosaic_dem_2tiles

3.)  เขียน GRADS control file เพื่อบอกให้โปรแกรมทราบ ตั้งชื่อว่า "ctl_read_merge_dem2tiles.ctl"
dset DEM_2tiles_merged.dat
title gtopo30 globaldem map 0.008333333333333 deg
undef -9999.0
XDEF 9600  LINEAR   60    0.00833333333333
YDEF 6000  LINEAR  -10    0.00833333333333
zdef 1 levels 0 1
tdef 1 linear 22Z24aug2010 1hr
vars 1
dem     0 99 elevation [m]
endvars
               
4)เปิด GRADS และเปิด control file ข้อสามมา
open ctl_read_merge_dem2tiles.ctl
set gxout shaded
d dem
cbarn

จะได้ผลลัพธ์ดังภาพ