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