All:
I figured out a way to do this interpolation, and Dennis Shea thought
I should share it on the message board. Improvements (or if you see
a bug) are welcome:
The problem was to interpolate from a cartesian regular grid to the
conic plane that a radar scans, but keeping a regular grid spacing
between points.
So the grid is actually equally spaced in x/y, but the height z of
the (x,y) is a function of distance from the radar and the elevation
angle.
I figure this out from the NCL examples page on interpolating to
isentropic coordinates. See http://www.ncl.ucar.edu/Applications/
Scripts/isent_1.ncl
Naturally, this was AFTER I sent the ncl-talk mail.
So the code looks something like this - I add in a lot of extraneous
coordinate information for debugging, and hopefully clarity.
;=======================================================================
===============
; this is the block of code is where you would read in your data, and
coordinate information
dbz = new((/nz,ny,nx/),float)
dbz = scalar field(:,:,:) ; some scalar field like
reflectivity, or u-wind component
xaxis = new((/nx/),float) ; array of x-coordinate values
yaxis = new((/ny/),float) ; array of y-coordinate values
zaxis = new((/nz/),float) ; array of z-coordinate values
xaxis = x_coord(:) ; axis data
yaxis = y_coord(:) ; axis data
zaxis = z_coord(:) ; axis data
; where is the radar relative to your coordinates and what is the
elevation angle you want to interpolate to and what is the height of
the radar relative to your z-coord?
radar_x_pos = -4500.
radar_y_pos = 4500.
radar_elev = 2.5 ; in degrees
radar_hgt = 30.0 ; meters, the physical
height of the radar dish
rcd = 45. / atan(1.0)
el = radar_elev / rcd
; Create 3D coordinate arrays, add in coordinate information for the
interpolation (int2p) program
xg = conform(dbz,xaxis,3)
yg = conform(dbz,yaxis,2)
zg = conform(dbz,zaxis,1)
zg!0 = "z"
zg!1 = "y"
zg!2 = "x"
zg&x = xaxis
zg&y = yaxis
zg&z = zaxis
; Create coordinates relative to radar position
xg = xg - radar_x_pos
yg = yg - radar_y_pos
; Compute height of radar surface - put into 3D array for int2p
z2d = radar_hgt + (xg(0,:,:)^2 + yg(0,:,:)^2)^0.5 * tan
(el)
z3d = new((/ny,nx,1/), float)
z3d(:,:,0) = z2d(0,:,:)
z3d!0 = "y"
z3d!1 = "x"
z3d!2 = "z"
z3d&x = xaxis
z3d&y = yaxis
; Call int2p, note that we reorder the axes because int2p wants the
interpolating coordinate to be the rightmost dimension - we already
created z3d with that coordinate order
dbz_radar_plane = int2p(zg(y|:,x|:,z|:), dbz(y|:,x|:,z|:), z3d,
0) ; output array has dimensions (/ny,nx,1/)
; For radial velocity, you need to interpolate the u,v,w fields to
the plane, then compute the azimuthal angles...
; You should also probably adjust the vertical velocity to add in the
contamination from the falling precip terminal, but we will ignore
that here
az = atan2(yg(0,:,:),xg(0,:,:)) ;
output array has dimensions (/ny,nx/)
radial_velocity = u_radar_plane(:,:,0)*cos(az) + v_radar_plane
(:,:,0)*sin(az) + w_radar_plane(:,:,0)*sin(el) ; output array has
dimensions (/ny,nx,1/)
;=======================================================================
================================================
Hope this is helpful.
regards,
Lou
------------------------------------------------------------------------
---- | Dr. Louis J. Wicker | Research Scientist, National Severe Storms Lab | 1313 Halley Cir, Norman, OK 73069 | E-mail: Louis.Wicker@noaa.gov | HTTP: www.nssl.noaa.gov/~lwicker | Phone: (405) 366-0416 | Fax: (405) 366-0472 | | | Two wrongs don't make a right, but three left's do... ------------------------------------------------------------------------ ---- | | "The contents of this message are mine personally and | do not reflect any position of the Government or NOAA." | ------------------------------------------------------------------------ ---- > >> >> All: >> >> I am looking for some hints on how to solve the following problem >> without using 2 for loops in NCL. >> >> The problem is I am trying to interpolate 3D gridded data fld to a >> radar conical surface. I simply want the xy grid to remain the same, >> but >> to interpolate to that 2D plane given by some elevation angle and >> distance from the radar. >> >> So I have a scalar 3D field like reflectivity >> >> dbz(z,y,x) >> >> which has xg, yg, and zg grid dimensions (each are 1D). >> >> and I have a 2D field of heights created from knowing where each xy >> grid pt is and computing the distance from the radar so that >> >> zhgt(y,x) = distance_from_radar * tan(elevation_angle). >> >> >> So at each x,y point, I need to search for the correct height point >> near zhgt(y,x) and then do a simple linear interpolation. >> >> The dumb way would be (forgive the bad syntax) >> >> for j = 0,ny-1 >> for i = 0,nx-1 >> >> k = (search for vertical grid pt right below zhgt(j,i)) >> >> uzinterp = coeff * dbz(k,j,i) + (1-coefff)*dbz(k+1,j,i) >> >> end for >> end for >> >> But I will bet for a 200x200 (or larger grid) this will be slow! >> >> I am pretty sure that I can easily create a 2D array of the vertical >> index "k" >> >> k_2d = ind(zg .le. zhgt) ???? >> >> but how can I select the 2D set of points in the dbz(z,y,x) that are >> given by k_2d? is there some fancy notation, or do I have to convert >> the >> dbz array to 1D, select the points out in some rather nasty fashion, >> and then convert back to 2D? >> >> Any help would be appreciated, as I am trying to get a nice figure >> done for the meso/radar conference next week.. >> >> Cheers! >> >> Lou Wicker >> >> >> _______________________________________________ >> ncl-talk mailing list >> ncl-talk@ucar.edu >> http://mailman.ucar.edu/mailman/listinfo/ncl-talk >> > >
_______________________________________________ ncl-talk mailing list ncl-talk@ucar.edu http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2b29 : Wed Oct 26 2005 - 18:25:11 MDT