2014.07.31
ขั้นตอน-การเฉลี่ยแบบจำลองความสูงด้วย Kernel
หรือ Moving
window
ผมต้องการจะเฉลี่ยค่าแบบจำลองความสูง GTOPO30 1ระวาง ด้วยพื้นที่โดยรอบของจุดที่พิจารณา ในที่นี้ใช้หลักการ Moving window แล้วทดลองทำแต่ละระยะทางจากจุดความสูงที่กำลังพิจารณาที่
10,25,50 เพื่อทำการเปรียบเทียบ และนำไปใช้ในงานอื่นต่อไป
1.
ใช้ Fortran
ในการเฉลี่ยตามโค้ดด้านล่าง
program
avg_gtopo30
!----------------------------------------------------------
! 2014.07.30
! read gtopo 1 tile after
swap to little endian
! apply average in horizontal length at
considering grid
! 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,k
!-average dem in horizonal length
integer
:: si,ei,sj,ej
integer
:: m,n
real,parameter
:: rgd=10.0 !-prefered
radius of grid [pixel]
real,dimension(ncol,nrow)
:: demsm !-smooth
dem output
real
:: tdem
real,parameter
:: min_elv=-1.0E+8 !- limit of
minimum elevation
real,parameter
:: max_elv= 1.0E+8 !- limit of
maximum elevation
real,parameter
:: na_elv=-9999
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)//'E100N40.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
print*,'smooth...'
do
i=1,ncol
do j=1,nrow
!-start and
end pixel of defined radius
si=i-rgd
ei=i+rgd
sj=j-rgd
ej=j+rgd
!print*,i,j,si,ei,sj,ej
!-unwanted
area
if ((i<=rgd) .or. ((ncol-i)<=rgd))
then
demsm(i,j)=demn(i,j)
cycle
endif
!-unwanted
area
if ((j<=rgd) .or. ((nrow-j)<=rgd))
then
demsm(i,j)=demn(i,j)
cycle
endif
!-loop for
average dem
tdem=0.0;k=0
do m=si,ei
do n=sj,ej
!-check
outlier in dem
if
((demn(m,n)>max_elv).or.(demn(m,n)
!print*,demn(m,n),'over!!'
cycle
endif
!print*,i,j,m,n,demn(m,n)
tdem=tdem+demn(m,n)
k=k+1
enddo
enddo
demsm(i,j)=tdem/k !-divided by total
number of pixel
!print*,demsm(i,j),demn(i,j),tdem,k
!print*,'----------------------------------------------------------'
enddo
enddo
open(unit=20,file='DEM_tile_E100N40_SM10.dat',form='unformatted',access='direct',recl=row*col*sizeofreal)
write(20,rec=1)demsm
close(20)
print*,'done...'
endprogram
avg_gtopo30
2.
สร้าง Control
file ของ Grads ชื่อ ctl_read_dem1tile_sm10.ctl
dset
DEM_tile_E100N40_SM10.dat
title
GTOPO30 map 0.008333333333333 deg
undef
-9999.0
**options little_endian template
**options little_endian
**options yrev 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
3.
แสดงผลใน Grads
open
ctl_read_dem1tile_sm10.ctl
set
gxout shaded
d
dem
cbarn
4.
ตัวอย่างประกอบด้านล่าง
ไม่มีความคิดเห็น:
แสดงความคิดเห็น