Fwd: [ncl-talk] Is there a way to do this without for loops? Solution!

From: Louis Wicker (Louis.Wicker AT noaa.gov)
Date: Tue Oct 25 2005 - 09:45:23 MDT

• Next message: Mike Notaro: "num"

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

radar_elev = 2.5 ; in degrees
radar_hgt = 30.0 ; meters, the physical

rcd = 45. / atan(1.0)

; 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

; 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

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/)

(:,:,0)*sin(az) + w_radar_plane(:,:,0)*sin(el) ; output array has
dimensions (/ny,nx,1/)
;=======================================================================
================================================

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
>>
>> 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