[Thread Prev][Thread Next][Index]

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,

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/
***************************************
[Thread Prev][Thread Next][Index]
Contact Us
Dept of Commerce / NOAA / OAR / PMEL / TMAP

Privacy Policy | Disclaimer | Accessibility Statement