2014.08.06
Squall line เป็นแนวฝนแบบ Convective
ที่จัดเรียงตัวกันในมาตราส่วน Meso scale คิดง่ายคือ
ประมาณหลักร้อยกิโลเมตรนั่นคือ
พื้นที่เรดาร์ตรวจอากาศภาคพื้นดินสามารถตรวจจับรูปแบบการเคลื่อนที่ของฝนนี้ได้
ในแถบอินโดจีนมีโปรเฟสเซอร์ Prof.Takehiko Satomura ได้ทำการอธิบายและเสนอกลไกการกำเนิดรูปแบบฝนชนิดไว้แล้ว
ในปี2000 น่าเสียดายที่ท่านได้จากโลกนี้ไปแล้วเมื่อเดือนมีนาคม 2014 ท่านเป็นอาจารย์ของผมเอง
หากแต่ไอเดียท่านยังอยู่กับพวกเรา
แนวคิดภาพตัดขวางแนวดิ่ง RHI
1.
อ่านเรดาร์แบบ Volume ไฟล์ และเลือกมุมอซิมุทที่ต้องการ Plot เลือกความสูง
ระยะทางของเรดาร์บีม
2.
คำนวนความสูงและระยะทาง ของแต่ละจุดที่เรดาห์ได้ทำการบันทึกข้อมูล
3.
เรียกซับโปรแกรมเพื่อประมาณค่าจุดเรดาร์ในแต่ละมุมยก ด้วย IDW
ตัวอย่างฝนชนิด Squalline
program
create_RHI
!-------------------------------------------------------------
! 2014.08.06 NATAPON MAHAVIK
! to create data for plotting Range-Height-Indicator [RHI]
! from extracted radar of
Phetchabun by Gematronik library
! data stor at
/home/mnattapon/RAID3storage4/research_icp/
!
analysis_phb_alg/radar_ans/cal_RHI_PPI_CAPPI_rad/cal_RHI
! data of selected azimuth will be
interpolate by IDW
!-------------------------------------------------------------
implicit none
integer,parameter ::
sizeofreal=4 !-record size
integer :: i,j,m,n
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,parameter :: k2m=1000.0 !- km to meters
real,parameter :: NA=-9999.0 !- no radar data
!-radar variables
integer :: rn !-record number
integer :: maxelv !-max number of
elevations
real :: elv_rad ! elevation
angle in radian
character*10 :: path
real :: rgate ! real(igate) = gate "km"
real :: rbin ! real(irbin) = rbin "km"
integer :: ielv, iazm, ibin
integer :: maxazmt !-Max number of
rays as temp var
integer:: maxbin ! Max number of
range
integer,parameter ::
maxazmf=1000 !-Max number of
rays as my asssumption!!!
real, dimension(:,:,:), allocatable
:: rdata !-Radar data
real, dimension(:,:),
allocatable :: elv !-Elevation angles
real, dimension(:,:),
allocatable :: azmth !-Angles [ in mathematic in degree]
real,dimension(:,:), allocatable
:: azmth_rad !-Angles [in azimuth in
radian]
integer,dimension(:),allocatable
:: maxazm !-Max number of rays
depending on sweeps
real,dimension(:,:,:),allocatable
::beam_dist_hz !-beam distance in
horizontal plain in m.
real,dimension(:,:,:),allocatable
::beam_height !-beam height by reinhart
p62 in m.
integer,parameter::re=6374000. !-earth radius in m.
real,parameter ::rr=4./3.*re !-reinhart p62 as R' in
m.
real,parameter :: rad_alt=97. !-radar altitude above MSL
from calculation DEM.(GG~73)
real :: h_rd=30. !-height of
radar antenna by estimate TMD using 30m
integer :: az_int !-round up real azimuth to interger [degree]
real,parameter :: sel_az=225.0 !-selected azimuth degree
real,parameter ::
min_dbz=10.0,max_dbz=100.0 !-max and
min dbz
real,parameter :: max_dist=200.0 !-max distance for RHI in km
real,parameter :: max_high=15.0 !-max hight for RHI in km
!-output var for interpolation
real, dimension(:,:), allocatable
:: dbz !-dbz at elevation(i) & bin number(i) [dBZ]
real, dimension(:,:), allocatable
:: height !-height from ellipsoid
at elevation(i) & bin number(i) [km]
real, dimension(:,:), allocatable
:: dist_x !-distance from radar
site at elevation(i) & bin number(i)[km]
real, parameter :: res=1.0 !- unit in km
integer,parameter :: nx=max_dist/res !- number of x pixels
integer,parameter ::
ny=max_high/res !- number of
y pixels
real, dimension(nx,ny)::
outp_intp !- interpolated
output
real, parameter :: radi=2.0 !-radius of interpolation
in idw [km]
!-read radar binary
!-cut
rn=1
open(10,file=trim(path)//'PHB1008102100.dat',form='unformatted',access='direct',recl=sizeofreal)
read(10,rec=rn) maxelv
!-cut
!-initialize output
dbz=NA;height=NA;dist_x=NA
do ielv = 1, maxelv
maxazmt=maxazm(ielv)
do
iazm = 1, maxazmt
read(10,rec=rn) azmth(iazm,ielv)
rn = rn + 1
!--- convert azmuth to be mathematic angle
!-cut
read(10,rec=rn) elv(iazm,ielv)
rn = rn + 1
do ibin = 1, maxbin
read(10,rec=rn) rdata(ibin,iazm,ielv)
rn = rn + 1
!-cal beam position in horizontal plain
!-cut
!-cal beam height + rad altitude in meter.
!-cut
endif
end do
end do
end do
close(10)
!-call sub IDW interpolate_idw(dbz(height,dist))
print*,'interpolate dbz to idw of
range-height radar ...'
call interpolate_idw(dbz,height,dist_x,maxelv,maxbin,res,nx,ny,radi,max_dbz,outp_intp)
print*,'infos for GRADS ;)
xdef:',nx,'ydef: ',ny,'resolution: ',res
!-write interpolated radar RHI to binary for checking in grads
open
(20,file='rs_interpolated_RHI_radar.bin',form='unformatted',access='direct',
recl=nx*ny*sizeofreal)
write(20,rec=1)outp_intp
close(20)
open
(30,file='rs_interpolated_RHI_radar_ncl.bin',form='unformatted')
write(30)outp_intp
close(30)
deallocate(beam_dist_hz,beam_height);deallocate(rdata,elv,azmth,azmth_rad)
deallocate(maxazm);deallocate(dbz,height,dist_x)
end
program create_RHI
Subroutine
interpolate_idw(dbz,height,dist_x,maxelv,maxbin,res,nx,ny,radi,max_dbz,outp_intp)
implicit none
!-external
integer,intent(in) ::
maxelv,maxbin
real,dimension(maxelv,maxbin),intent(in)
:: dbz,height,dist_x
real,intent(in) :: max_dbz
real,intent(in) :: res !-defined resolution of output grid
[km]
integer :: nx,ny
real,dimension(nx,ny),intent(out)
:: outp_intp
real :: radi !-defined radius for interpolation [km]
!-internal variables
real :: dist !-distance between output grid
and input grid
integer :: sx,ex,sy,ey !-start and end grid of considering box
respected to radius
integer :: i,j,m,n
real,parameter :: NA=-9999.0 !-not
available data
real,parameter :: p=-2 !-power for idw
integer :: x,y !-coordinate of output grid in
input data grid
real :: inv_d !-inverse distance
real :: nw !-normalized weight
real,dimension(nx,ny) ::
sum_invd !-sum inverse distance
real,dimension(nx,ny) ::
sum_nw !-sum normalized weight
!-loop elvation and bin of input data to be fast calculate
!outp_intp=NA
!-1st loop for find summation of inverse distance in output grid
do i=1,maxelv
do j=1,maxbin
if (dbz(i,j)==NA)cycle
!-find box of considering interpolation respected to defined radius
sx=((int(dist_x(i,j))-radi)/res)-res; ex=((int(dist_x(i,j))+radi)/res)+res
sy=((int(height(i,j))-radi)/res)-res; ey=((int(height(i,j))+radi)/res)+res
!-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
!-loop inside box for calculate summation of inverse distance
do m=sx,ex
do n=sy,ey
!-convert coordinate considering
grid to be same input coordinates
x=(m*res)-res+1; y=(n*res)-res+1
!-calculate euclidean distance
dist=sqrt(abs(x-dist_x(i,j))+abs(y-height(i,j)))
!print*,x,y,dist
!-cal inverse distance and sum all
inverse distance
inv_d=dist**p
sum_invd(m,n)=sum_invd(m,n)+inv_d
enddo
enddo
enddo
enddo
!-2nd loop apply normalized weight to dbz and sum all results of each
ouput grid
do i=1,maxelv
do j=1,maxbin
if (dbz(i,j)==NA)cycle
!-find box of considering interpolation
respected to defined radius
sx=((int(dist_x(i,j))-radi)/res)-res; ex=((int(dist_x(i,j))+radi)/res)+res
sy=((int(height(i,j))-radi)/res)-res; ey=((int(height(i,j))+radi)/res)+res
!-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
!-loop inside box for calculate weight, then apply it
do m=sx,ex
do n=sy,ey
!-convert coordinate considering grid to be same input
coordinates
x=(m*res)-res+1; y=(n*res)-res+1
!-calculate euclidean distance
dist=sqrt(abs(x-dist_x(i,j))+abs(y-height(i,j)))
!-cal inverse distance, normalize and sum all weighted
dbz
inv_d=dist**p
nw=inv_d/sum_invd(m,n)
sum_nw(m,n)=sum_nw(m,n)+nw !-it must be 1.0 in every
grid
outp_intp(m,n)=outp_intp(m,n)+(nw*dbz(i,j))
enddo
enddo
enddo
enddo
!-replace values for grads or NCL plot
where (outp_intp>max_dbz)
outp_intp=0
endwhere
End subroutine interpolate_idw
5.
ตัวอย่างโค้ดของ NCL .ในการ Plot ภาพ RHI โปรแกรมนี้มีอะไรน่าตื่นเต้นเยอะนะครับ
ช่วยในงานวิจัยได้ไม่น้อยเลยครับ
;**********************************************************************
; color_1.ncl
; NATTAPON MAHAVIK
; Concepts illustrated:
;
- Drawing color filled contours using the default color map
; how to read and write data from
f90 to ncl
;
https://www.ncl.ucar.edu/Document/Functions/Built-in/fbinread.shtml
;
;**********************************************************************
load
"$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load
"$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;************************************************
begin
;************************************************
; read in binary
;************************************************
fili="rs_interpolated_RHI_radar_ncl.bin"
recl=200*15
;rad =
fbindirread("rs_interpolated_RHI_radar_ncl.bin",rec,(/15,200/),"float")
r=fbinread (fili, recl , "float")
rad=onedtond(r,(/15,200/))
rad@_FillValue = -9999 ;missing values
;************************************************
; create plot
;************************************************
wks = gsn_open_wks("ps","ps_RHI_radar_PHB") ; open a ps file
gsn_define_colormap(wks,"prcp_1") ; choose color map
res =
True ; plot mods desired
res@vpWidthF = 0.8 ; set width and height
res@vpHeightF = 0.5
res@cnFillOn =
True ; turn on color fill
res@gsnMaximize = True
res@gsnDraw =
False ; don't draw
yet
res@gsnFrame =
False ; don't
advance frame yet
res@cnInfoLabelOn =
False ; not show infomation height
res@cnLineLabelsOn =
False ; not show label
of contour line
res@cnConstFLabelOn = False
res@cnLineThicknessF = 0 ; change line thickness
res@cnLinesOn = False
res@tiMainString =
"Range Height Indicator of Squallline in Thailand"
res@tiXAxisString =
"Distance [km]"
res@tiYAxisString =
"Height[km]"
res@cnLevelSelectionMode = "ManualLevels" ; manually set the contour levels with
the following 3 resources
res@cnMinLevelValF = 10.0 ; set the minimum contour
level
res@cnMaxLevelValF = 60.0 ; set the maximum contour
level
res@cnLevelSpacingF = 5.0 ; contour interval
res@lbOrientation = "vertical" ; vertical label bar
res@lbLabelStride = 2.0
res@gsnSpreadColors=True
plot = gsn_csm_contour(wks,rad,res)
draw(plot) ; note
we are drawing the first one!
frame(wks)
end
สิ่งที่ต้องทำต่อไปคือ แก้โปรแกรมให้ยืดหยุ่น ใช้กับข้อมูลหลายช่วงเวลา นำวิธีการประมาณค่าแบบอื่นๆมาใช้
ไม่มีความคิดเห็น:
แสดงความคิดเห็น