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/
***************************************
use test.nc define grid/like=test grd ! list edges list yboxlo[g=grd],yboxhi[g=grd] ! calculate weights w1 and w2 let rad = 3.14159265/180 let w1 = (yboxhi[i=1,g=grd]-yboxlo[i=1,g=grd])*cos(y[i=1,g=grd]*rad) let w2 = sin(yboxhi[i=1,g=grd]*rad) - sin(yboxlo[i=1,g=grd]*rad) ! list normalized weights list w1/w1[j=@sum],w2/w2[j=@sum] ! list weights calculated by Ferret at j=1,2,3,4,5,18 let all = if test gt 0 then 1 let total_area = all[x=@din,y=@din] let cell1 = if test eq 1 then 1 let cell2 = if test eq 2 then 1 let cell3 = if test eq 3 then 1 let cell4 = if test eq 4 then 1 let cell5 = if test eq 5 then 1 let cell18 = if test eq 18 then 1 list cell1[x=@din,y=@din]/total_area list cell2[x=@din,y=@din]/total_area list cell3[x=@din,y=@din]/total_area list cell4[x=@din,y=@din]/total_area list cell5[x=@din,y=@din]/total_area list cell18[x=@din,y=@din]/total_area
Attachment:
test.nc
Description: Cdf file