Re: [ferret_users] spatial regridding

[Ansley Manke <ansley.b.manke@xxxxxxxx>]

Hi Emilie,
Ferret does do the average you want. As you say, it's a matter of lining 
up the input grid and the output grid. The sum you give,
( E + (D+B+F+H)/2 + (A+C+G+I)/4 ) / 4

The sum that Ferret computes is the sum of the variable value at the 
coordinate locations, weighted by the portion of the grid box 
surrounding that coordinate which lies within the output grid cell.  If 
the new grid has its coordinate locations at every second coordinate 
point from the first grid, then the output grid box contains all of grid 
box E, half of grid boxes B, D, H, and F, and a quarter of grid boxes A, 
C, G, and I.

A  B  C
D  E  F
G  H  I

There is an FAQ on a related topic, 
http://ferret.pmel.noaa.gov/Ferret/FAQ/analysis/regridding_with_ave.html 
which you might look at. The ideas there could be extended to more than 
one dimension. Here is a simple example, which will also illustrate what 
is done.

Let's define a simple grid and define a variable on it

yes? def axis/x=0.5:6.5:1 x7
yes? def axis/y=0.5:6.5:1 y7
yes? let var = cos(x[gx=x7]) + sin(y[gy=y7])

! List commands will show full single-precision values
yes? set list /prec=6 
yes? list var
            VARIABLE : COS(X[GX=X7]) + SIN(Y[GY=Y7])
            SUBSET   : 7 by 7 points (X-Y)
             0.5      1.5      2.5      3.5      4.5      5.5      6.5  
             1        2        3        4        5        6        7
0.5 / 1:  1.35701  0.55016 -0.32172 -0.45703  0.26863  1.18810  1.45601
1.5 / 2:  1.87508  1.06823  0.19635  0.06104  0.78670  1.70616  1.97408
2.5 / 3:  1.47605  0.66921 -0.20267 -0.33798  0.38768  1.30714  1.57506
3.5 / 4:  0.52680 -0.28005 -1.15193 -1.28724 -0.56158  0.35789  0.62580
4.5 / 5: -0.09995 -0.90679 -1.77867 -1.91399 -1.18833 -0.26886 -0.00094
5.5 / 6:  0.17204 -0.63480 -1.50668 -1.64200 -0.91634  0.00313  0.27105
6.5 / 7:  1.09270  0.28586 -0.58602 -0.72134  0.00432  0.92379  1.19171

! Now define a new grid, with the coordinates at 1.5, 3.5, and 5.5 in 
both directions

def axis/x=1.5:5.5:2 x3
def axis/y=1.5:5.5:2 y3
def grid/x=x3/y=y3 xy3

yes?  list var[g=xy3@ave]
            VARIABLE : COS(X[GX=X7]) + SIN(Y[GY=Y7])
                       regrid: XY3 on X@AVE, on Y@AVE
            SUBSET   : 3 by 3 points (X-Y)
             1.5      3.5      5.5  
             1        2        3
1.5 / 1:  0.82270  0.04701  1.31400
3.5 / 2: -0.21568 -0.99137  0.27563
5.5 / 3: -0.48889 -1.26459  0.00241

! Now, compare the result of the regridding at output point (1.5, 1.5) 
with the result of your ABCDEF formula,

yes? list (1.06823 + (0.55016 + 1.87508 + 0.66921 + 0.19635)/2 \
       +   (1.35701+1.47605-0.20267-0.32172)/4 ) /4

            VARIABLE : (1.06823 + (0.55016+1.87508+0.66921+0.19635)/2 + 
(1.35701+1.47605-0.20267-0.32172)/4 ) /4
         0.822699

Which matches the output regridding at (1.5,1.5)

To define the output grid so that the coordinates line up, you may want 
to see the coordinate values of you input grid,

yes? list x[gx=my_var]
yes? list y[gy=my_var]

So, if you want to define an output grid whose grid cells are twice as 
large in each direction with every second coordinate from the input 
grid, you will get the result from your equation. Ferret always does 
this kind of weighted average, however, even if the cells are not lined 
up as I did here, or if the coordinates of one or both grids are 
irregularly spaced.

If the units of the axes are longitude and latitude, Ferret applies a 
cosine(latitude) correction to the computed averages.

Emilie Vanvyve wrote:
> Hi,
>
> Could anybody help me to understand how exactly Ferret computes a 
> spatial average?
>
> I have a regular horizontal grid, its resolution is 40 km and I am 
> trying to reduce the resolution to 80, 160 and 320 km. The initial x 
> axis ranges from -2880 to 2440 km and the y axis from -2560 to 2760 
> km. The coordinates of each grid box correspond to the centre of the 
> box: for example, the grid box centred in (-2880, -2560) covers the 
> area x=-2900:-2860 and y=-2580:-2540, i.e. a square of 1600 km2.
>
> The computation works fine and the resolution is reduced nicely—I get 
> the same using @SBX:2 or @SBX:4 or @SBX:8, with the difference that in 
> the latter case the plot resolution is still 40 km.
> The problems arises when I try to check how Ferret did the 
> computation, I get totally lost and I assume there is probably a trick 
> I don't know.
> I know how work the time axis and that the time step must correspond 
> to the middle of the average period in order to take into account the 
> correct time steps. Is there something similar with the spatial axes?
>
> Here is what I do to go from the 40-km resolution to the 80-km 
> resolution (initially in a script) :
>
> yes? define axis/x=-2840:2360:80/unit=km xax
> yes? define axis/y=-2520:2680:80/unit=km yax
> yes? define grid/x=xax/y=yax grd
> yes? let var80km = var40km[g=grd@ave]
>
> The new x and y axis are slightly shorter so that no missing points 
> are found at the boundaries.
> What I want Ferret to do is the following. Let's zoom on 9 grid boxes:
> A  B  C
> D  E  F
> G  H  I
>
> In order to get the 80-km resolution, E should be replaced by: ( E + 
> (D+B+F+H)/2 + (A+C+G+I)/4 ) / 4, meaning that D, B, F and H 
> contributes for half of their surface, A, C, G and I for a quarter. 
> When I check the values of var80km, it is not what I get. Why?
>
> Thanks very much for any idea on this!
>
> Emilie V.
--

~>   ~>   ~>   ~>   ~>   ~>   ~>   ~>   ~>   ~>   ~> 
Ansley Manke, NOAA/PMEL   ansley.b.manke@noaa.gov
7600 Sand Point Way NE, Seattle WA
Phone 206-526-6246,  FAX 206-526-6744