การจัดแต่งโครงสร้างโมเดลความสูง 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
ได้ภาพดังข้างล่างนี้
ไม่มีความคิดเห็น:
แสดงความคิดเห็น