**Previous message:**Louis Wicker: "Is there a way to do this without for loops?"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ]

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

**Next message:**Mike Notaro: "num"**Previous message:**Louis Wicker: "Is there a way to do this without for loops?"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ]

*
This archive was generated by hypermail 2b29
: Wed Oct 26 2005 - 18:25:11 MDT
*