Re: linint2_points for station data

From: Dennis Shea (shea AT cgd.ucar.edu)
Date: Thu Nov 18 2004 - 15:39:54 MST

  • Next message: Gary Strand: ""debug mode" for NCL?"

    >
    >Hello, a couple of weeks ago I sought help on extracting
    >point data from a grid, given a lat/lon that didn't necessarily
    >correspond to the grid points. I was led to linint2_points() and
    >have made "some" progress, but of course, I had to complicate
    >the initial requirements a little bit :)
    >
    >Given a regular 3x3 grid varying from 45N to 47N and 112W to 114W,
    >and a temperature field, I can easily extract/interpolate point
    >data as follows:
    >
    >;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    >begin
    >
    > temperature = (/ (/100, 100, 50/), \
    > (/100, 100, 100/), \
    > (/100, 100, 100/) /)
    >
    > lat = (/45, 46, 47/)
    > lon = (/-114, -113, -112/)
    >
    > ; Here's the lat/lon I want station data for
    > y0 = 45.5
    > x0 = -112.5
    >
    > station_temperature = linint2_points(lon, lat, temperature, \
    > False, x0, y0, 0)
    >
    > print(station_temperature)
    >
    >end
    >;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    >
    >However, what if my grid is NOT so "regular?" Suppose, as is the
    >case with some of my WRF output, a "row" of points might start at
    >45N and end at 45.2N? Does linint2_points() allow for this? In
    >the following, I've naively tried by placing all of my lat/lon grid
    >points in arrays and passing those to linint2_points(), as follows
    >
    >;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    >begin
    >
    > temperature = (/ (/100, 100, 100/), \
    > (/100, 50, 100/), \
    > (/100, 100, 100/) /)
    >
    > lat_points = (/ (/45, 45.1, 45.2 /), \
    > (/46, 46.1, 46.2 /), \
    > (/47, 47.1, 47.2 /) /)
    >
    > lon_points = (/ (/-114, -113, -112/), \
    > (/-114.1, -113.1, -112.1/), \
    > (/-114.2, -113.2, -112.2/) /)
    >
    > ; Here's the lat/lon of the location I want station data for
    > y0 = 45.5
    > x0 = -112.5
    >
    > station_temperature = linint2_points(lon_points, lat_points, \
    > temperature, False, x0, y0, 0)
    >
    > print(station_temperature)
    >
    >end
    >;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    >
    >
    >but I get an error message saying:
    >
    >fatal:linint2_points: If xi is not one-dimensional, then it must have
    >one less dimension than fi
    >
    >(Note that above, "xi" corresponds to "lon_points" and "fi" corresponds
    >to "temperature").
    >
    >Indeed, the documentation confirms that xi must have one less dimension
    >than fi, but
    >I don't understand why this would be the case (which suggests I don't
    >really
    >understand the purpose of multi-dimension arguments for xi and yi). Is
    >there any place
    >I can go for a little more explanation? A Google search on
    >"linint2_points" doesn't
    >turn up much :)

    ============
    Hello

    linint2_points was developed to interpolate from
    a "regular" grid to arbitrary locations within the regular grid.
    By "regular", I mean a grid where any point on the grid
    can be accessed via one-dimensional (1d) coordinate variables.
    To be considered a "coordinate variabe" it must be a 1D
    array of monotonically {in/de}creasing values. They need not
    be equally spaced.

    The WRF grid is curvilinear. It requires 2D lat/lon arrays to
    describe the grid. It is not so easy to develop a fast, general
    purpose code to interpolate from any arbitrary grid to a point.
    Suffice it to say, some grids are *very* unusual.

    I asked the MMM/WRF 'user support' person if they had
    an interpolation function. Here is the response:

    > We probably do not have the exact thing you
    > ask for. We have some routines that converts
    > x, y to lat, long, and vice versa. But the x,y
    > here refers to the coordinate on projected maps.
    > If that would be useful to you, I will pass them to you.

    ---
    NCL's Natgrid interpolation library may work.
    I have never used Natgrid. Also, I have heard it can
    be slow. The following is untested 
    

    http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/nnptinit.html http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/nnpnt.html

    Let's say lat2d(lat,lon) and lon2d(lat,lon) are the grid coordinates. for the variable x(ntim,lat,lon). Let latSta(N), lonSta(N), xSta(ntim,N) be the station observations

    lat1d = ndtooned(lat2d) lon1d = ndtooned(lon2d) do nt=0,ntim-1 nnpntinit(lon1d, lat1d, ndtooned(x(nt,:,:)) ) xSta = nnpnts(lonSta,latSta) ; xSta(N) : end do Good Luck D

    _______________________________________________ ncl-talk mailing list ncl-talk AT ucar.edu http://mailman.ucar.edu/mailman/listinfo/ncl-talk



    This archive was generated by hypermail 2b29 : Thu Nov 18 2004 - 17:10:20 MST