# Re: [ferret_users] Inaccurate area weights calculation

```Hi all,
```
Yes, we should be able to change this for future releases. It certainly helps to have these examples with large grid cells, where the discrepancy is largest. Thank you for the report, Marco.
```
Ansley

Andrew Wittenberg wrote:
```
```Hi Marco,

```
I think the problem may be that Ferret uses the following formula to compute 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 is not the case -- but if we redefine your cell centers to be midway between your 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 be independent of where the cell center is located within a grid cell. [Ansley, I've reopened Trac #1348 with a request to replace the first formula above with the second. Presumably this wouldn't be too difficult?]
```
Andrew

On Mon, 17 Aug 2009, Marco Steinacher wrote:

```
```Dear Andrew,

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

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

```
```
```