Re: [ferret_users] Inaccurate area weights calculation

To: Marco Steinacher <steinacher@xxxxxxxxxxxxxxxx>
Subject: Re: [ferret_users] Inaccurate area weights calculation
From: Andrew Wittenberg <Andrew.Wittenberg@xxxxxxxx>
Date: Mon, 17 Aug 2009 17:40:31 -0400 (EDT)
Cc: oar.pmel.ferret_users@xxxxxxxx
In-reply-to: <4A899715.6000603@xxxxxxxxxxxxxxxx>
References: <4A89608D.8070408@xxxxxxxxxxxxxxxx> <Pine.LNX.4.64.0908171154590.28226@atw> <4A899715.6000603@xxxxxxxxxxxxxxxx>
Reply-to: Andrew Wittenberg <Andrew.Wittenberg@xxxxxxxx>
Sender: owner-ferret_users@xxxxxxxx

Hi Marco,

I think the problem may be that Ferret uses the following formula tocompute the cell areas (yc = cell center, yh = cell height):


   cos(yc) * 2*abs(sin(yh/2))

which is only identical to

   abs(sin(yhi)-sin(ylo))

when yc lies exactly midway between ylo and yhi. In your example, this isnot the case -- but if we redefine your cell centers to be midway betweenyour cell edges, with


def ax/y/bounds/unit=degrees_n `test,r=yaxis` = \
   (yboxhi[g=test]+yboxlo[g=test])/2, \
   yboxlo[g=test], yboxhi[g=test]

then we get the answers you expected for the weightings.

Thank you for pointing this out; the area calculation really ought to beindependent of where the cell center is located within a grid cell.[Ansley, I've reopened Trac #1348 with a request to replace the firstformula above with the second. Presumably this wouldn't be toodifficult?]


Andrew


On Mon, 17 Aug 2009, Marco Steinacher wrote:

Dear Andrew,

Thanks for your reply. I'm using the following version:

FERRET v6.2
Linux(g77) 2.4.21-32 - 05/19/09

From your message I infer that this version of Ferret should use
weighting w2 or something similar. I have attached two files to test
this for my specific grid. The (shorted) output is:

$ ferret -script test.jnl
Column  1: YBOXLO is YBOXLO (axis LAT_T)
Column  2: YBOXHI is YBOXHI (axis LAT_T)
          YBOXLO  YBOXHI
76.5S /  1: -90.00 -70.81
66.4S /  2: -70.81 -62.73
59.4S /  3: -62.73 -56.44
53.7S /  4: -56.44 -51.06
48.6S /  5: -51.06 -46.24
[...]
1.6S  / 18:  -3.18   0.00
[...]

Column  1: EX#1 is W1/W1[J=@SUM]
Column  2: EX#2 is W2/W2[J=@SUM]
              EX#1    EX#2
76.5S /  1:  0.03826  0.02778
66.4S /  2:  0.02750  0.02778
59.4S /  3:  0.02725  0.02778
53.7S /  4:  0.02718  0.02778
48.6S /  5:  0.02716  0.02778
[...]
1.6S  / 18:  0.02712  0.02778

VARIABLE : CELL1[X=@DIN,Y=@DIN]/TOTAL_AREA
0.03811
VARIABLE : CELL2[X=@DIN,Y=@DIN]/TOTAL_AREA
0.02750
VARIABLE : CELL3[X=@DIN,Y=@DIN]/TOTAL_AREA
0.02725
VARIABLE : CELL4[X=@DIN,Y=@DIN]/TOTAL_AREA
0.02719
VARIABLE : CELL5[X=@DIN,Y=@DIN]/TOTAL_AREA
0.02716
VARIABLE : CELL18[X=@DIN,Y=@DIN]/TOTAL_AREA
0.02713


I think this shows that Ferret is using something similar to weighting
w1 (although not exactly the same) instead of the correct weighting w2.
Or am I missing something?

Many thanks,
Marco


Andrew Wittenberg wrote:

Dear Marco,

What version of Ferret are you running?  The area weighting should be
correct for Ferret versions 6.00 and later (i.e. computed as your w2).
[Developers, refer to Trac ticket #1348, fixed April 2007.]

Andrew


On Mon, 17 Aug 2009, Marco Steinacher wrote:

Dear Ferret Users,

I'm struggling with a problem related to the inaccurate calculation of
weights for area-weighted operations like @ave and @din in Ferret when
using a very coarse grid.

The grid has 36x36 equal-area cells and it is regular spaced in
longitude (10deg). In the NetCDF file the coordinates of the center
(lat_t) and of the edges (lat_u using the 'edges' attribute) is given
(see below).

When looking only at the grid cells at one longitude (eg. i=1) the
weight of each grid cell should be 1/36 = 0.028. However, the weights
calculated by Ferret differ largely. Obviously Ferret calculates the
weights as

 w1 = (lat_u[j+1]-lat_u[j])*dlon*cos(lat_t[j])

instead of

 w2 = ( sin(lat_u[j+1]) - sin(lat_u[j]) )*dlon

(Here, dlon is the spacing of the grid cells in longitude and, of
course, everything is converted to radians for the calculation). For j=1
we get the following (normalized) weights

 w1 = 0.038
 w2 = 0.028.

w2 is correct but the weight calculated by Ferret is more than 30%
higher!

I'm wondering if anybody has experienced similar problems and if there
is any workaround to force Ferret to use the correct area-weights for
@ave etc. Any comments and suggestions on this issue are very welcome!

Many thanks,
Marco


PS: I know that I can use var[x=@sum,y=@sum]/var[x=@ngd,y=@ngd] to get
the arithmetic average which is equal to the weighted average in this
specific case. Nevertheless, I would prefer to be able to use @ave, @din
etc. as well.

--------------------------------------------------------

double lat_t(lat_t) ;
               lat_t:long_name = "Latitude (T grid)" ;
               lat_t:units = "degrees_north" ;
               lat_t:minimum = -76.4637972626189 ;
               lat_t:maximum = 76.4637972626188 ;
               lat_t:edges = "lat_u" ;

double lat_u(lat_u) ;
               lat_u:long_name = "Latitude (U grid)" ;
               lat_u:units = "degrees_north" ;
               lat_u:minimum = -90. ;
               lat_u:maximum = 89.9999991462264 ;

lat_t = -76.4637972626189, -66.4435356908988, -59.4415682140651,
       -53.6639424853861, -48.5903778907291, -43.9829631303587,
       -39.7090165530684, -35.6853347126521, -31.8554308240259,
       -28.1786427405299, -24.6243183521641, -21.1684488458325,
       -17.7915905730076, -14.4775121859299, -11.2122714176497,
        -7.98355614555541, -4.78019184719916, -1.59175417658591,
         1.5917541765859,   4.78019184719916,  7.98355614555541,
        11.2122714176497,  14.4775121859299,  17.7915905730076,
        21.1684488458324,  24.6243183521641,  28.1786427405299,
        31.8554308240259,  35.6853347126521,  39.7090165530684,
        43.9829631303587,  48.5903778907291,  53.6639424853861,
        59.4415682140651,  66.4435356908988,  76.4637972626188 ;

lat_u = -90,               -70.8118635462791, -62.7339555492672,
       -56.4426902380793, -51.0575587310186, -46.2382573073202,
       -41.8103148957786, -37.6698869643296, -33.7489885958886,
       -30,               -26.387799961243,  -22.8853804761586,
       -19.4712206344907, -16.1276202131608, -12.8395884069042,
        -9.5940682268606,  -6.37937020844281, -3.18473853672041,
        -3.18055468146e-15, 3.18473853672041,  6.3793702084428,
         9.59406822686046, 12.8395884069041,  16.1276202131608,
        19.4712206344907,  22.8853804761586,  26.387799961243,
        30,                33.7489885958886,  37.6698869643296,
        41.8103148957786,  46.2382573073202,  51.0575587310186,
        56.4426902380793,  62.7339555492672,  70.8118635462791,
        89.9999991462264 ;

--
***************************************
Marco Steinacher

Climate and Environmental Physics
Physics Institute, University of Bern
Sidlerstrasse 5, CH-3012 Bern

Phone ++41 (0)31 631 34 02
Fax   ++41 (0)31 631 87 42
steinacher@xxxxxxxxxxxxxxxx
http://www.climate.unibe.ch/
***************************************



--
***************************************
Marco Steinacher

Climate and Environmental Physics
Physics Institute, University of Bern
Sidlerstrasse 5, CH-3012 Bern

Phone ++41 (0)31 631 34 02
Fax   ++41 (0)31 631 87 42
steinacher@xxxxxxxxxxxxxxxx
http://www.climate.unibe.ch/
***************************************

Follow-Ups:
- Re: [ferret_users] Inaccurate area weights calculation
  - From: Ansley Manke

References:
- [ferret_users] Inaccurate area weights calculation
  - From: Marco Steinacher
- Re: [ferret_users] Inaccurate area weights calculation
  - From: Andrew Wittenberg
- Re: [ferret_users] Inaccurate area weights calculation
  - From: Marco Steinacher

Previous by thread: Re: [ferret_users] Inaccurate area weights calculation
Next by thread: Re: [ferret_users] Inaccurate area weights calculation

[Thread Prev][Thread Next][Index]

Dept of Commerce / NOAA / OAR / PMEL / TMAP