ขั้นตอนการประมาณค่าข้อมูลลมECWMF
ด้วยวิธี IDW
ข้อมูลลมจากหน่วยงาน ECWMF มีความละเอียดที่
1.5 องศาต่อหนึ่งจุดภาพ
หากต้องการได้ข้อมูลที่มีความละเอียดกว่านี้ต้องทำการประมาณค่าเองซึ่งมีหลายวิธีการตั้งแต่ง่ายไปจนถึงซับซ้อน
แล้วหากพื้นที่ศึกษามีขนาดใหญ่แล้วช่วงเวลาที่ใช้ในการศึกษาปรากฎการณ์มีช่วงเวลายาวซึ่งใช้ในงานภูมิอากาศวิทยา
ทำให้ต้องเขียนโปรแกรมเพื่อให้ได้ข้อมูลที่ต้องการ
ผมได้ทำการทดลองเขียนโปรแกรมฟอร์แทรน(ท่านผู้อ่านอยากจะใช้อะไรก็ได้นะครับ C, MATLAB, Python, IDL)เพื่อประมาณค่าให้ได้ข้อมูลละเอียดขึ้นเป็น 0.5 องศาด้วยวิธี IDW อาศัยข้อมูลหลักจากมุมทั้งสี่ที่ใกล้เคียงกับจุดที่กำลังประมาณค่าอยู่ตามภาพนะครับ หมายความว่าหนึ่งบล๊อกจะมีค่าที่ต้องประมาณอยู่9ค่า โปรแกรมอาจจะยาวหน่อย คงต้องพัฒนาให้มันสั้นลงกว่านี้ครับ แล้วลองการประมาณการแบบอื่นๆ หากพบข้อผิดพลาดรบกวนบอกกันด้วยนะครับ
1.
ดาวน์โหลดข้อมูลจากเวปไซต์ ECMWF เอาลมผิวพื้น
10 เมตรทั้งสองคอมโพเนนท์คือ U,V
2.
ใช้โปรแกรม WGRIBเพื่อแปลงเป็นข้อมูล
ไบนารีไฟล์
3.
พัฒนาโปรแกรมดังโค้ดข้างล่าง
program simple_weighted_interpolation
!-------------------------------------------------------------------------------------------
! 2014.08.02 NATTAPON MAHAVIK
! interpolation
of coarse grid wind ECMWF data [1.5degree] to finer grid
! here we use
inverse weighted distance of
! those grid
to coarse grid
! concept of interpolation
IDW:
! https://www.ems-i.com/smshelp/Data_Module/Interpolation/Inverse_Distance_Weighted.htm
!---------------------------------------------------------------------------------------------------
implicit none
!-variable for
ECMWF
integer,parameter
:: ncol=29,nrow=35,ntim=4 !-x and y
dimension know it when we download
real,parameter
:: slon=99.0,slat=-10.5
!-lower left of
wind data dowloaded from ECMWF
degree
real,parameter
:: rese=1.5 !-coarse grid resolution of
ECMWF degree
real,dimension(ncol,nrow)
:: u,v,un,vn !-wind orographic
and moisture convergence flux
integer ::
i,j,m,n
integer,
parameter :: shft_xu=1, shft_xv=3 !-shifting in
x-axis of u,v by comparing to origial data [pixels]
!-variable for
output grid
real,parameter
:: lat1=slat,lon1=slon,lat2=slat+(nrow*rese)-rese,lon2=slon+(ncol*rese)-rese !-boundary in geographic coordinate of coarse data
real,parameter
:: res=0.5
!-prefered fine grid resolution in degree ctl
file GRADs
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),parameter
:: wu_int=0.,wv_int=0. !-array of
result here is wind
real,dimension(4):: wu !-wind U of coarse grid at any corners
[m/s]
real,dimension(4):: wv !-wind V of coarse grid at any corners
[m/s]
integer,parameter
:: sizeofreal=4
real :: ln1,ln2,lt1,lt2 !-geographic corrdinate at
any 4 grid corners
print*,'number
of column and row of new fine grid for ctl GRADS: ',nx,ny
print*,lat1,lat2,lon1,lon2
print*,'resolution:
',res,' degree'
print*,'interpolation
by IDW...'
!--read
wind data got from WGRIB
open (10,file='../wind_surface_2010081012utc_EPR_ICP.bin',form='unformatted',access='direct',recl=ncol*nrow*sizeofreal)
read(10,rec=1)u
read(10,rec=2)v
close(10)
!-inverse
wind data and shift
do i=1,ncol
do j=1,nrow
m=ncol-i+1 !-col inverse
n=nrow-j+1 !-row inverse
un(i-shft_xu,n)=u(i,j)
vn(i-shft_xv,n)=v(i,j)
enddo
enddo
!-call
interpolation fucntion IDW of each grid corrner
do i=1,ncol
do j=1,nrow
!--find geographic
coordinates of coarse grid at any corners
ln1=slon+(i*rese)-rese
ln2=slon+((i+1)*rese)-rese
lt1=slat+(j*rese)-rese
lt2=slat+((j+1)*rese)-rese
!print*,lt1,lt2,ln1,ln2
!-find u and v as respected
to any 4 corner grids
wu=(/un(i,j),un(i+1,j),un(i+1,j+1),un(i,j+1)/)
wv=(/vn(i,j),vn(i+1,j),vn(i+1,j+1),vn(i,j+1)/)
!--call interpolation
function IDW
call intp_by_dist(slon,slat,lt1,lt2,ln1,ln2,res,nx,ny,wu,wv,wu_int,wv_int)
enddo
enddo
!--write wind
U,V to binary for checking in grads
open (20,file='wind_int_wind_surface_2010081012_ICP.bin',form='unformatted',access='direct',
recl=nx*ny*sizeofreal)
write(20,rec=1)wu_int
write(20,rec=2)wv_int
close(20)
endprogram simple_weighted_interpolation
subroutine intp_by_dist(slon,slat,lat1,lat2,lon1,lon2,res,nx,ny,wu,wv,wu_int,wv_int)
!-----------------------------------------------------------------
! calculate weigth by distance and apply to varaibles [wind]
!-----------------------------------------------------------------
implicit none
real,intent(in)
:: slon,slat
real,intent(in)
:: lat1,lon1,lat2,lon2,res
!-coordinates of any 4 corner grids
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,m,n
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 :: inv_d1,inv_d2,inv_d3,inv_d4 !- inverse distance
real ::
sum_inv_dist
real,parameter
:: p=-2 !- power of
inverse distance
integer ::
sx,sy,ex,ey !- start and end of lon and
lat axis for 4 corner
grids[grid]
real,parameter
:: NA=-9999.0
real,parameter
:: wind_err1=1000.0,wind_err2=-1000.0
!-Function
calculating distance on sphere
real ::
haversine !-Function name for calculating
distance on sphere
!print*,'interpolaing
by inverse distance weighted...'
!---------------------------------------------------------------------
!-box location of 4 corner grids in new fine grid [important step!!]
sx=int((lon1-slon)/res)+1
ex=int((lon2-slon)/res)+1
if ((lat1<=0) .and. (lat2<=0)) then
sy=int((abs(slat)-abs(lat1))/res)+1
ey=int((abs(slat)-abs(lat2))/res)+1
endif
if ((lat1<=0) .and. (lat2>=0)) then
!sy=int(abs(slat-abs(lat1))/res)+1
!ey=int((slat+lat2)/res)+1
sy=int(abs(abs(slat)-abs(lat1))/res)+1
ey=int((abs(slat)+lat2)/res)+1
endif
if ((lat1>=0) .and. (lat2>=0)) then
sy=int((abs(lat1)+abs(slat))/res)+1
ey=int((abs(lat2)+abs(slat))/res)+1
endif
!---------------------------------------------------------------------
!-here to avoid
error outbound otherwise segmentation fault
if (sx<=0) sx=1
if (ex>nx)
ex=nx
if (sy<=0) sy=1
if (ey>ny)
ey=ny
!-interpolate
only inside any 4 grid
corners
do m=sx,ex
do n=sy,ey
!-cal considering geocoordinate
lon=slon+((m*res)-res)
lat=slat+((n*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)
!-inverse distance
inv_d1=(d1**p);inv_d2=(d2**p);inv_d3=(d3**p);inv_d4=(d4**p)
if (d1==0.0) inv_d1=0.0; if (d2==0.0) inv_d2=0.0 ;if (d3==0.0) inv_d3=0.0;if (d4==0.0) inv_d4=0.0
sum_inv_dist=inv_d1+inv_d2+inv_d3+inv_d4
nw1=inv_d1/sum_inv_dist
nw2=inv_d2/sum_inv_dist
nw3=inv_d3/sum_inv_dist
nw4=inv_d4/sum_inv_dist
!-set highest weight at
corrner grid location
if (d1==0.0) then
nw1=1.0;nw2=0.0;nw3=0.0;nw4=0.0
endif
if (d2==0.0) then
nw2=1.0;nw1=0.0;nw3=0.0;nw4=0.0
endif
if (d3==0.0) then
nw1=0.0;nw2=0.0;nw3=1.0;nw4=0.0
endif
if (d4==0.0) then
nw1=0.0;nw2=0.0;nw3=0.0;nw4=1.0
endif
!-apply normalized weight to
each wind value
wu_int(m,n)=(wu(1)*nw1) + (wu(2)*nw2) + (wu(3)*nw3) + (wu(4)*nw4)
wv_int(m,n)=(wv(1)*nw1) + (wv(2)*nw2) + (wv(3)*nw3) + (wv(4)*nw4)
!-to avoid interpolated
error caused by input data
if ((wu_int(m,n)>wind_err1).or.(wv_int(m,n)>wind_err1).or.(wu_int(m,n)2).or.(wv_int(m,n)2)) then
wu_int(m,n)=NA
wv_int(m,n)=NA
endif
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
!-set criteria in latitude [important]
if ((dlat1<=0.0) .and. (dlat2<0.0)) dlat=abs(dlat1-dlat2)*d2r
if ((dlat1<=0.0) .and. (dlat2>0.0)) dlat=(abs(dlat1)+dlat2)*d2r
if ((dlat1>=0.0) .and. (dlat2>=0.0)) dlat=(abs(dlat1-dlat2))*d2r
if ((dlat1>=0.0) .and. (dlat2<=0.0)) dlat=(abs(dlat2)+dlat1)*d2r
dlon=abs(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
4.
ทำการตรวจสอบในGRADS แต่ต้องเขียน Control file
ก่อนดังข้างล่างนะครับ
dset
^wind_int_wind_surface_2010081012_ICP.bin
title
Interpolated sample wind map 0.5 deg
undef -9999.0
XDEF 85 LINEAR
99.0 0.5
YDEF 103 LINEAR
-10.5 0.5
zdef 1 levels 0 1
tdef 1 linear 12Z10aug2010 1hr
vars 2
u 0 99 zonal wind [m/s]
v 0 99 meridional wind [m/s]
จะได้ผลดังภาพข้างล่างเปรียบเทียบระหว่างก่อนทำการประมาณค่ากับหลังประมาณค่า หากจะให้ได้ข้อมูลที่ละเอียดกว่านี้ก็เพียงเปลี่ยนตรง res นะครับ
ไม่มีความคิดเห็น:
แสดงความคิดเห็น