2014.08.03
ขั้นตอนการทำการประมาณค่าน้ำฝนด้วยวิธี IDW จากข้อมูลตัวอย่าง
การประมาณค่าข้อมูลเชิงพื้นที่ใดๆที่ได้มาจากภาคสนามนั้นเป็นขั้นตอนสำคัญต่อการวิเคราะห์ข้อมูลเชิงพื้นที่ จากข้อมูลที่กระจัดกระจายในพื้นที่ นำมาผ่านวิธีการประมาณค่าให้เป็นข้อมูลเชิงพื้นผิวหรือข้อมูลเชิงกริดเพื่อให้ง่ายต่อการวิเคราะห์ วิธีการที่นำมาใช้ในการประมาณมีตั้งแต่ง่ายไปจนถึงซับซ้อน ไม่ว่าจะเป็น linear interpolation, Thiessen, Fixed-radius, Natural neighbor, Trend surface interpolation, IDW, Spline, Kriging เป็นต้น ซึ่งในแต่ละประเภทก็จะมีแยกย่อยไปเยอะแยะ แล้วแต่ความต้องการในงานนั้นๆ
ในส่วนของผมนั้น ผมต้องการทำการประมาณค่าน้ำฝนรายชั่วโมงหลายช่วงฤดูกาล
หากใช้โปรแกรมช่วยพวก GIS, R, Python หรือโปรแกรมเชิงพื้นที่อื่นๆก็อาจจะต้องใช้เทคนิคเข้าช่วยเพราะข้อมูลเยอะ
ในที่นี้ผมจึงอยากทดลองเขียนโปรแกรมด้วยฟอร์แทรนเพื่อการประมวลผลข้อมูลน้ำฝนตัวอย่าง
หากมีข้อมูลเยอะสามารถใช้ bash script
เข้าช่วยในการประมวลผลแล้วเพียงนั่งรอทำงานอื่นไป โปรแกรมที่นำมาแสดงนี้ใช้วิธีการ
IDW โดยสามารถกำหนดระยะรัศมีของการประมาณค่า และค่ายกกำลังของส่วนกลับระยะทาง
ดูภาพประกอบนะครับ โปรแกรมจะคำนวนค่าน้ำหนักตามส่วนกลับของระยะทางในรัศมีที่ได้ระบุไว้
ลองเอาไปเล่นดูนะครับ หากมีข้อบกพร่องประการใดแจ้งให้ทราบได้เลยนะครับ
จะได้ปรับปรุงต่อไป
1.
โปรแกรมฟอร์แทรน IDW (Inverse Distance Weighted Interpolation)
program
simple_weighted_interpolation
!------------------------------------------------------------------------
! 2014.08.03 NATTAPON MAHAVIK
! interpolation of scattered data (i.e.rain) to fine grid
! suppose we have set of rain at stations that are located over observation area
! we want to interpolate it simply to fine grid.
! here we use inverse weighted distance of
! those scatteredd points to fine grid at prefered radius distance
! concept of interpolation IDW:
!------------------------------------------------------------------------
implicit none
real,parameter :: lat1=10.0,lon1=100.0,lat2=40.0,lon2=140.0 !-geographic coordinate of boundary interpolated data
real,parameter :: res=1.0 !-grid resolution in degree
real,parameter :: radi=20.0 !-prefered radius of interpolation degree
integer,parameter :: nx=nint((lon2-lon1+res)/res),ny=nint((lat2-lat1+res)/res) !-number of grid in x,y direction
real*4,dimension(nx,ny) :: outp_int !-array of result here is rainfall rate
!-rain sample data
integer,parameter :: nums=15 !-number of rain stations
real,dimension(nums),parameter:: rain=(/10.0,5.0,1.0,20.0,30.0,10.0,5.0,1.0,20.0,30.0,10.0,5.0,1.0,20.0,30.0/) !-rain rate [mm]
real,dimension(nums),parameter:: lonr=(/105.5,110.0,115.0,127.0,135.0,105.5,110.0,115.0,127.0,135.0,105.5,110.0,115.0,127.0,135.0/) !-longitude of rain stations [degree]
real,dimension(nums),parameter:: latr=(/15.0,13.0,17.0,19.0,14.0,25.0,23.0,27.0,29.0,24.0,35.0,33.0,37.0,39.0,34.0/) !-latitude of rain stations
integer,parameter :: sizeofreal=4
print*,'number of column and row: ',nx,ny
print*,'resolution: ',res,' degree'
!------------------------------------------------------------------------
! 2014.08.03 NATTAPON MAHAVIK
! interpolation of scattered data (i.e.rain) to fine grid
! suppose we have set of rain at stations that are located over observation area
! we want to interpolate it simply to fine grid.
! here we use inverse weighted distance of
! those scatteredd points to fine grid at prefered radius distance
! concept of interpolation IDW:
!------------------------------------------------------------------------
implicit none
real,parameter :: lat1=10.0,lon1=100.0,lat2=40.0,lon2=140.0 !-geographic coordinate of boundary interpolated data
real,parameter :: res=1.0 !-grid resolution in degree
real,parameter :: radi=20.0 !-prefered radius of interpolation degree
integer,parameter :: nx=nint((lon2-lon1+res)/res),ny=nint((lat2-lat1+res)/res) !-number of grid in x,y direction
real*4,dimension(nx,ny) :: outp_int !-array of result here is rainfall rate
!-rain sample data
integer,parameter :: nums=15 !-number of rain stations
real,dimension(nums),parameter:: rain=(/10.0,5.0,1.0,20.0,30.0,10.0,5.0,1.0,20.0,30.0,10.0,5.0,1.0,20.0,30.0/) !-rain rate [mm]
real,dimension(nums),parameter:: lonr=(/105.5,110.0,115.0,127.0,135.0,105.5,110.0,115.0,127.0,135.0,105.5,110.0,115.0,127.0,135.0/) !-longitude of rain stations [degree]
real,dimension(nums),parameter:: latr=(/15.0,13.0,17.0,19.0,14.0,25.0,23.0,27.0,29.0,24.0,35.0,33.0,37.0,39.0,34.0/) !-latitude of rain stations
integer,parameter :: sizeofreal=4
print*,'number of column and row: ',nx,ny
print*,'resolution: ',res,' degree'
!--call interpolation function
IDW
call intp_by_dist(lat1,lat2,lon1,lon2,res,nx,ny,radi,nums,rain,latr,lonr,outp_int)
!--write interpolated rain to binary for checking in grads
open (10,file='rain_int_sample.bin',access='direct', recl=nx*ny*sizeofreal)
write(10,rec=1)outp_int
close(10)
endprogram simple_weighted_interpolation
subroutine intp_by_dist(lat1,lat2,lon1,lon2,res,nx,ny,radi,nums,rain,latr,lonr,outp_int)
!-------------------------------------------------------
! calculate weigth by distance and apply to varaibles
implicit none
real,intent(in) :: lat1,lon1,lat2,lon2,res,radi
integer,intent(in) :: nx,ny,nums
real,dimension(nums),intent(in) :: rain
real,dimension(nums),intent(in) :: latr
real,dimension(nums),intent(in) :: lonr
real,dimension(nx,ny),intent(out):: outp_int
real,dimension(nx,ny):: sum_inv_dist
real,dimension(nx,ny):: sum_weight_rain
integer :: sx,ex,sy,ey ! start and end grid within defined radius of interpolation
real :: lon,lat ! geo coordinate of considering grid
integer :: i,j,m,n
real :: dist !-distance between considering rain sta. and fine grid in [m]
real,parameter :: p=-2 !- power of inverse distance
!-Function calculating distance on sphere
real :: haversine !-Function name for calculating distance on sphere
print*,'interpolaing by inverse distance weighted...'
!-loop each rain station
do i=1,nums
sx=nint(((lonr(i)-radi)-lon1)/res)
ex=nint(((lonr(i)+radi)-lon1)/res)
sy=nint(((latr(i)-radi)-lat1)/res)
ey=nint(((latr(i)+radi)-lat1)/res)
!-loop within defined radius in output grid
do m=sx,ex
do n=sy,ey
!-geocoordinate of considering grid
lon=lon1+(m*res)
lat=lat1+(n*res)
!-to avoid out of boundary
if ((m<1 m="" n="" or.="">nx).or.(n>ny)) cycle
!-cal distance of considering rain and output grid
dist=haversine(lat,lon,latr(i),lonr(i))
if (dist==0.0) then
sum_inv_dist(m,n)=1.0
sum_weight_rain(m,n)=rain(i)
cycle
else
sum_inv_dist(m,n)=sum_inv_dist(m,n)+(dist**p)
sum_weight_rain(m,n)=sum_weight_rain(m,n)+((dist**p)*rain(i))
end if
end do
end do
end do
!-evaluate interpolated rain by IDW
outp_int=sum_weight_rain/sum_inv_dist
end subroutine intp_by_dist
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 :: d2r=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)*d2r
dlon=(dlon2-dlon1)*d2r
lat1=dlat1*d2r
lat2=dlat2*d2r
!-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
สร้างคอนโทรลไฟล์ของโปรแกรม GRADS ตามข้างล่างแล้วเปิดตรวจดูใน GRADS
call intp_by_dist(lat1,lat2,lon1,lon2,res,nx,ny,radi,nums,rain,latr,lonr,outp_int)
!--write interpolated rain to binary for checking in grads
open (10,file='rain_int_sample.bin',access='direct', recl=nx*ny*sizeofreal)
write(10,rec=1)outp_int
close(10)
endprogram simple_weighted_interpolation
subroutine intp_by_dist(lat1,lat2,lon1,lon2,res,nx,ny,radi,nums,rain,latr,lonr,outp_int)
!-------------------------------------------------------
! calculate weigth by distance and apply to varaibles
implicit none
real,intent(in) :: lat1,lon1,lat2,lon2,res,radi
integer,intent(in) :: nx,ny,nums
real,dimension(nums),intent(in) :: rain
real,dimension(nums),intent(in) :: latr
real,dimension(nums),intent(in) :: lonr
real,dimension(nx,ny),intent(out):: outp_int
real,dimension(nx,ny):: sum_inv_dist
real,dimension(nx,ny):: sum_weight_rain
integer :: sx,ex,sy,ey ! start and end grid within defined radius of interpolation
real :: lon,lat ! geo coordinate of considering grid
integer :: i,j,m,n
real :: dist !-distance between considering rain sta. and fine grid in [m]
real,parameter :: p=-2 !- power of inverse distance
!-Function calculating distance on sphere
real :: haversine !-Function name for calculating distance on sphere
print*,'interpolaing by inverse distance weighted...'
!-loop each rain station
do i=1,nums
sx=nint(((lonr(i)-radi)-lon1)/res)
ex=nint(((lonr(i)+radi)-lon1)/res)
sy=nint(((latr(i)-radi)-lat1)/res)
ey=nint(((latr(i)+radi)-lat1)/res)
!-loop within defined radius in output grid
do m=sx,ex
do n=sy,ey
!-geocoordinate of considering grid
lon=lon1+(m*res)
lat=lat1+(n*res)
!-to avoid out of boundary
if ((m<1 m="" n="" or.="">nx).or.(n>ny)) cycle
!-cal distance of considering rain and output grid
dist=haversine(lat,lon,latr(i),lonr(i))
if (dist==0.0) then
sum_inv_dist(m,n)=1.0
sum_weight_rain(m,n)=rain(i)
cycle
else
sum_inv_dist(m,n)=sum_inv_dist(m,n)+(dist**p)
sum_weight_rain(m,n)=sum_weight_rain(m,n)+((dist**p)*rain(i))
end if
end do
end do
end do
!-evaluate interpolated rain by IDW
outp_int=sum_weight_rain/sum_inv_dist
end subroutine intp_by_dist
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 :: d2r=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)*d2r
dlon=(dlon2-dlon1)*d2r
lat1=dlat1*d2r
lat2=dlat2*d2r
!-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
สร้างคอนโทรลไฟล์ของโปรแกรม GRADS ตามข้างล่างแล้วเปิดตรวจดูใน GRADS
dset
^rain_int_sample.bin
title Interpolated sample wind map 1.0 deg
undef -9999.0
XDEF 41 LINEAR 100 1.000
YDEF 31 LINEAR 10 1.000
zdef 1 levels 0 1
tdef 1 linear 22Z24aug2010 1hr
vars 1
r 0 99 rainfall rate interpolated [mm]
endvars
title Interpolated sample wind map 1.0 deg
undef -9999.0
XDEF 41 LINEAR 100 1.000
YDEF 31 LINEAR 10 1.000
zdef 1 levels 0 1
tdef 1 linear 22Z24aug2010 1hr
vars 1
r 0 99 rainfall rate interpolated [mm]
endvars
ผลลัพธ์ที่ได้ตามภาพข้างล่างเป็นไปตามที่คาดไว้
ถือว่าโปรแกรมน่าจะไม่ผิดพลาดอะไรนะครับ หากต่อไปคงต้องพิจารณาจำนวนจุดที่ใช้ในการประมาณค่าแทนรัศมี หรือทดลองการประมาณค่าแบบที่ซับซ้อนมากขึ้น
ไม่มีความคิดเห็น:
แสดงความคิดเห็น