วันพฤหัสบดีที่ 31 กรกฎาคม พ.ศ. 2557

ขั้นตอน-การเฉลี่ยแบบจำลองความสูงด้วย Kernel หรือ Moving window



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.                  ตัวอย่างประกอบด้านล่าง
 

ไม่มีความคิดเห็น:

แสดงความคิดเห็น