เป็นวิธีการต่อโมเดลความสูงของผมเอง โดยผมต้องการต่อเพื่อเอาไปประมวลผลความสูงกับลมพื้นผิวในขั้นต่อไปเพื่อหาฝนที่เกี่ยวเนื่องกับภูเขา โปรแกรมอาจจะไม่สวยเพราะเขียนแบบลูกทุ่ง แต่ก็ใช้งานได้ ลองเอาไปทำกันดูครับ
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
จะได้ผลลัพธ์ดังภาพ
ไม่มีความคิดเห็น:
แสดงความคิดเห็น