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"


    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

    I figure this out from the NCL examples page on interpolating to
    isentropic coordinates. See http://www.ncl.ucar.edu/Applications/
    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
         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.




    |  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