2014.08.02
ขั้นตอน-การประมาณค่าข้อมูลลมรายละเอียดสูงจากข้อมูลสมมติ
ผมต้องการที่จะประมาณค่าข้อมูลลมจากข้อมูลมูลที่หยาบโดยสมมติว่ามีข้อมูลอยู่
4 กริดห่างกัน แล้วต้องการประมาณค่าข้อมูลลมในลงไปในกริดย่อยๆนั้นที่รายละเอียด
1.0 องศาต่อหนึ่งกริด ในการทดลองนี้ผมใช้วิธีการ Inverse Distance Weighted Interpolation (IDW)ในการประมาณค่าข้อมูลลมที่หยาบนั้นๆ แล้วได้ตรวจสอบใน
GRADS พบว่าผลที่ได้เป็นไปตามที่คาดไว้
แสดงว่าโปรแกรมไม่น่าจะผิดนะ ^^’ ขั้นตอนตามนี้นะครับ
1.
โค้ดฟอร์แทรนสมมติลม U และ V
ตามมุมทั้งสี่แล้ว สร้างกริดย่อย 1 องศาเข้าไป แล้ววนลูปคำนวนระยะทางด้วย Haversine จากนั้นเอาระยะทางที่ได้มาคำนวนค่าน้ำหนัก ด้วย IDW
program simple_weighted_interpolation
!-------------------------------------------------------------------------------------------
! 2014.08.02 NATTAPON MAHVIK
! interpolation
of coarse grid data to fine grid
! suppose we
have 4 wind data at
the corner grids
! we want to
interpolate it simply to fine grid.
! here we use
inverse weighted distance of
! those grid
to coarse grid
! concept of
interpolation IDW:
!-------------------------------------------------------------------------------------------
implicit none
real,parameter
:: lat1=10.0,lon1=100.0,lat2=40.0,lon2=140.0 !-geographic coordinate of coarse data
real,parameter
:: res=1.0 !-grid
resolution in 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)
:: wu_int,wv_int
!-array of result here is wind
real,dimension(4),parameter::
wu=(/10.0,20.0,5.0,1.0/) !-wind U of coarse grid [m]
real,dimension(4),parameter::
wv=(/10.0,3.0,8.0,10.0/) !-wind V of coarse grid [m]
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,wu,wv,wu_int,wv_int)
!--write wind
U,V to binary for checking in grads
open (10,file='wind_int_sample.bin',access='direct',
recl=nx*ny*sizeofreal)
write(10,rec=1)wu_int
write(10,rec=2)wv_int
close(10)
endprogram
simple_weighted_interpolation
subroutine intp_by_dist(lat1,lat2,lon1,lon2,res,nx,ny,wu,wv,wu_int,wv_int)
!-------------------------------------------------------
! calculate
weigth by distance and apply to varaibles
implicit none
real,intent(in)
:: lat1,lon1,lat2,lon2,res
integer,intent(in)
:: nx,ny
real,dimension(4),intent(in)
:: wu,wv
real,dimension(nx,ny),intent(out)::
wu_int,wv_int
real ::
lon,lat !
geo coordinate of considering grid
integer ::
i,j
real :: d1,d2,d3,d4 !-distance
between considering grid and coarse grid in [m]
real :: w1,w2,w3,w4 !- weight by
inverse distance
real :: nw1,nw2,nw3,nw4 !- normalized
weight
real :: totw !-
summation of weight before normalize
real ::
dall
!- total distance from considering grid
real ::
sum_inv_dist
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...'
do i=1,nx
do j=1,ny
!-cal considering
geocoordinate
lon=lon1+((i*res)-res)
lat=lat1+((j*res)-res)
!-cal dist function
d1=haversine(lat,lon,lat1,lon1)
d2=haversine(lat,lon,lat1,lon2)
d3=haversine(lat,lon,lat2,lon2)
d4=haversine(lat,lon,lat2,lon1)
sum_inv_dist=(d1**p)+(d2**p)+(d3**p)+(d4**p)
nw1=(d1**p)/sum_inv_dist
nw2=(d2**p)/sum_inv_dist
nw3=(d3**p)/sum_inv_dist
nw4=(d4**p)/sum_inv_dist
if (d1==0.0) nw1=1.0; if (d2==0.0) nw2=1.0 ;if (d3==0.0) nw3=1.0;if (d4==0.0) nw4=1.0
!-apply normalized
weight to each wind value
wu_int(i,j)=(wu(1)*nw1) + (wu(2)*nw2) + (wu(3)*nw3) + (wu(4)*nw4)
wv_int(i,j)=(wv(1)*nw1) + (wv(2)*nw2) + (wv(3)*nw3) + (wv(4)*nw4)
!print*,'-->',lon,lat,d1,d2,d3,d4,nw1,nw2,nw3,nw4,wu_int(i,j),wv_int(i,j)
enddo
enddo
end
subroutine intp_by_dist
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 :: 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
2.
สร้าง Control file ของ GRADS ดังนี้ “ctl_int_sample_wind.ctl”
dset
^wind_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 2
u 0 99 zonal wind [m/s]
v 0 99 meridional wind [m/s]
endvars
3.
ตรวจสอบใน GRADS
open ctl_int_sample_wind.ctl
set gxout
shaded
d u
d u;v
cbarn
ผลที่ได้ดังภาพข้างล่างเป็นไปตามที่คาดไว้
แสดงว่าโอเคครับ
ไม่มีความคิดเห็น:
แสดงความคิดเห็น