Geographie Freiheit: ขั้นตอน-การหาอัตราการเปลี่ยนแปลงความสูงของแบบจำลอง GTOPO30

2014.08.02

ขั้นตอนการหาอัตราการเปลี่ยนแปลงความสูงของแบบจำลอง GTOPO30

ผมต้องการหาอัตราการเปลี่ยนแปลงความสูง (Gradient)ของแบบจำลอง GTOPO30 ในทิศทางแกน X (Zonal)และ Y (Meridional) ซึ่งให้สอดคล้องกับข้อมูลลม เพื่อจะนำไปประมวลผลในขั้นต่อไป ขั้นตอนหลักๆคือ เฉลี่ย DEM หาGradientทั้งสองแกนในแต่ละค่าจุดภาพ ส่งออกผลเป็นไบนารี่เพื่อแสดงใน Grads มาเริ่มกันเลย ยาวหน่อยครับ แถมเรื่องการคำนวนระยะทางบนพื้นผิวโลกให้ด้วยโดยสูตร Haversine ในฟังค์ชั่นนะครับ เด็กที่ชอบเรื่องแผนที่ทราบไว้ก็จะเยี่ยมเลย

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

2. ใช้ฟอร์แทรนเฉลี่ย DEM จากนั้นหาGradientทั้งสองแกนในแต่ละค่าจุดภาพ และเก็บเป็นไบนารี่ตามข้างล่าง (ยาวมาก แบบลูกทุ่งเลย) อย่าลืมเปลี่ยนขนาดของ Moving window ที่จะใช้เฉลี่ย DEM ด้วย หากใช้เยอะจะเสียเวลามาก เช่นใช้ 10 pixel จะใช้เวลาประมาณ 5 นาทีแต่ละเครื่องไม่เท่ากัน แล้วแต่สมรรถนะ

program cal_dem_gradient

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

! 2014.08.01 NATTAPON MAHAVIK

! Smooth DEM, Find Gradient in zonal and Meridional, write to binanry

! Please change grd=10.0 as you like

! Haversine codes :: http://rosettacode.org/wiki/Haversine_formula#Fortran

! Haversine formula :: http://www.ig.utexas.edu/outreach/googleearth/latlong.html

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

implicit none

!-Variables for DEM

character*150 :: infile_path

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

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

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

integer*2,dimension(col,row) :: deml,demr !-dem left and right tiles

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

integer :: i,j,k,m,n

real,dimension(col,row) :: demz,demm !- gradient of DEM in zonal and meridional

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=-9999

real,parameter :: res=0.008333333333333 !- DEM resolution in degree

real :: lat1,lon1 !- geographic coordinate of considering pixel

real :: lat2,lon2 !- geographic coordinate of neighboring pixel

real,parameter :: ul_lon=100.0 !- geographic coordinate of ul tile

real,parameter :: ul_lat=40.0 !- geographic coordinate of ul tile

integer,parameter :: rgd=1 !- grid space [m]

real :: grdz,grdm !- gradient in zonal and meridional

real :: delv !- different of elevation in gradient

!-Variables for smooth DEM

integer :: si,ei,sj,ej

real,parameter :: grd=10.0 !-prefered radius of grid from considering pixel[pixel] change!!!

real,dimension(col,row) :: demsm !-smooth dem output

real :: tdem

!-declare variable

real :: dist !-distance in meters

real :: slope !-atan(rise/run) in %

!-Function calculating distance on sphere

real :: haversine !-Function name for calculating distance on sphere

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=col*row*sizeofint)

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

read(10,rec=1)demr

close(10)

print*,'reading...'

do i=1,col

do j=1,row

m=col-i+1 !-col inverse

n=row-j+1 !-row inverse

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

enddo

print*,'smooth DEM...'

do i=1,col

do j=1,row

!-start and end pixel of defined radius

si=i-grd; ei=i+grd

sj=j-grd; ej=j+grd

!print*,i,j,si,ei,sj,ej

!-unwanted area

if ((i<=grd) .or. ((col-i)<=grd)) then

demsm(i,j)=na

cycle

endif

!-unwanted area

if ((j<=grd) .or. ((row-j)<=grd)) then

demsm(i,j)=na

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)

cycle

endif

tdem=tdem+demn(m,n)

k=k+1

enddo

demsm(i,j)=tdem/k !-divided by total number of pixel

!print*,'----------------------------------------------------------'

enddo

print*,'calculating gradient of DEM ...'

do i=1,col

do j=1,row

!-unwanted area

if ((i<=rgd) .or. ((col-i)<=rgd)) then

demm(i,j)=na;demz(i,j)=na

cycle

endif

!-unwanted area

if ((j<=rgd) .or. ((row-j)<=rgd)) then

demm(i,j)=na;demz(i,j)=na

cycle

endif

!-calculate zonal gradient

lat1=ul_lat-(j*res)

lon1=ul_lon+((i+1)*res)

lat2=ul_lat-(j*res)

lon2=ul_lon+((i-1)*res)

dist=haversine(lat1,lon1,lat2,lon2) !-call function

delv=demsm(i+1,j)-demsm(i-1,j)

grdz=delv/dist

!-calculate meridional gradient

lat1=ul_lat-(j+1*res)

lon1=ul_lon+(i*res)

lat2=ul_lat-(j-1*res)

lon2=ul_lon+(i*res)

dist=haversine(lat1,lon1,lat2,lon2) !-call function

delv=demsm(i,j+1)-demsm(i,j-1)

grdm=delv/dist

!-assign gradient to zonal and meridional

demm(i,j)=grdm

demz(i,j)=grdz

enddo

!-write binary file of zonal and meridional gradient for display in grads

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

write(20,rec=1)demz

write(20,rec=2)demm

close(20)

print*,'done...'

endprogram cal_dem_gradient

Real function haversine(dlat1,dlon1,dlat2,dlon2)

implicit none

real,intent(in) :: dlat1,dlon1,dlat2,dlon2 !- lat lon in degree converted to radian

real :: haversine !- distance will be return by this var [m]

real :: lat1,lat2

real :: a,c

real, parameter :: pi = 4*atan(1.0) !- exploit intrinsic atan to generate pi

real,parameter :: rd=pi/180.0 !- radian conversion degree to radian

real :: dlat,dlon !- difference in lat and lon in radian

real,parameter :: re=6372.8*1000.0 !- earth radius in m

dlat=(dlat2-dlat1)*rd

dlon=(dlon2-dlon1)*rd

lat1=dlat1*rd

lat2=dlat2*rd

!-Haversine formula

a=(sin(dlat/2))**2+cos(lat1)*cos(lat2)*(sin(dlon/2)**2)

c=2*atan2(sqrt(a),sqrt(1-a))

haversine=re*c

end function haversine

3. สร้าง Control file ชื่อ “ctl_cal_dem_smooth_gradient_SM10.ctl”

dset ^DEM_gradient_tile_E100N40_SM10.dat

title Gradient calculated from GTOPO30 map 0.008333333333333 deg

undef -9999.0

XDEF 4800 LINEAR 100 0.00833333333333

YDEF 6000 LINEAR -10 0.00833333333333

zdef 1 levels 0 1

tdef 1 linear 22Z24aug2010 1hr

vars 2

grdz 0 99 zonal gradient [unitless]

grdm 0 99 meridional gradient [unitless]

endvars

4. เปิด GRADS แล้วลอง plot ภาพดูจะได้ดังภาพตัวอย่างข้างล่าง เป็นการเปรียบเทียบอัตราการเปลี่ยนแปลงความสูงแต่ละระดับ Moving window จะเห็นว่าที่ไต้หวันเมื่อใช้ 10 pixel จะเห็นการเปลี่ยนแปลงความสูงชัดเจนคือสีแดงมีการลาดขึ้นส่วนสีน้ำเงินมีการลาดลงของภูมิประเทศในแนวแกน X หรือ Zonal direction

แถมให้นิดนึง ผมเปลี่ยนค่าสีตามโค้ดข้างล่างนะครับ ลองทำเล่นๆกันดู++

Open ctl_cal_dem_smooth_gradient_SM10.ctl

d demz

set rgb 16 0 0 255

set rgb 17 55 55 255

set rgb 18 110 110 255

set rgb 19 165 165 255

set rgb 20 220 220 255

set rgb 21 255 220 220

set rgb 22 255 165 165

set rgb 23 255 110 110

set rgb 24 255 55 55

set rgb 25 255 0 0

set gxout shaded

set clevs -0.05 -0.04 -0.03 -0.02 -0.01 0.01 0.02 0.03 0.04 0.05

set ccols 16 17 18 19 20 1 21 22 23 24 25

cbarn

Geographie Freiheit

วันเสาร์ที่ 2 สิงหาคม พ.ศ. 2557

ขั้นตอน-การหาอัตราการเปลี่ยนแปลงความสูงของแบบจำลอง GTOPO30

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

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

เกี่ยวกับฉัน

คลังบทความของบล็อก