[Thread Prev][Thread Next][Index]

Re: [ferret_users] Re: conservation of surface values according to @DIN operator




Hi Francois,

The reason that your integration comes to 1.735 is that your grid GPC1 is dimensionless. The axes need to be defined with the correct units
in order to generate the correct areas.

set axis/units="degrees_E" lon
set axis/units="degrees_N" lat

will fix things.

Cheers,
Russ

On 28/10/16 01:12, Francois Delclaux wrote:
Hi Russ,

I tried what you proposed, but unfortunately, it din't work.

There is something I don't understand

- I have one mask file mask_15s.nc , 15" resolution = 0.0041667, with only 0 and 1
 GRID GAP1
 name       axis              # pts   start                end
 LONGITUDE LONGITUDE         1020mr   75.502E(75.502)      79.748E(79.748)
 LATITUDE  LATITUDE           480 r   41.419N              43.415N
 LEVEL     Z (level)            1 r-  0                    0
 TIME      T (time)             1 r   0                    0


- I have one  precipitation file p_1mn.nc, resolution is 1mn = 0.0166667
GRID GPC1
name axis #  pts start end
LON    X     206     r 75.917 79.333
LAT    Y      76     r 41.75 43
LEVEL1 Z       1     r  0   0
TIME1 TIME 13149 r 01-JAN-1979 12:00 31-DEC-2014 12:00

- first computation of mark area at 15"
list mask[d=1, x=@din, y=@din]
result is 1.581E+10 in m2 . OK, it is the real surface


- second computation of mask area with regriding to 1mn with the grid of precipitation file GPC1
let ma1 = mask[d=1, gxy = GPC1@NRST]
with NRST operator, I have always 0 or 1 in my map.

the the result of list ma1[x=@din, y=@din] is 1.735 !!!


- third computation of mask area with regriding to 1mn with explicit definition of grid, the same as precip file:
define axis/x=75.917E:79.333E:0.016666 ax
define axis/y=41.75N:43N:0.016666 ay
define grid/x=ax/y=ay ggxy

GRID GGXY
name axis #        pts    start    end
AX    LONGITUDE    206mr 75.917E 79.334E(79.334)
AY    LATITUDE     76 r    41.75N 43N

then regriding on ggxy :
let ma2 = mask[d=1,gxy=ggxy@NRST]
yes? list ma2[x=@din, y=@din] = 1.578E+10

this result is OK.

So the result  @DIN operator seems to depend on the way by which grid is defined, either from a file, either from direct definition.

Regards

François


----- Le 26 Oct 16, à 1:32, Russ Fiedler <russell.fiedler@xxxxxxxx> a écrit :
Hi Francois,
 
 It's not the integration that's causing the problem here but the
 bilinear interpolation from the fine grid to the coarse grid (GIN1).
 @LIN is the default interpolation and we don't want it here.
 
 You'll typically have many fine grid points within the bounds of a
 coarse cell. Now, if the four points contributing to the interpolation
 to the fine grid are all 1
 then, of course, you get 1 as the value on the coarse grid signifying
 that all points are good. BUT, if there are any masked points within the
 cell's boundaries then
 this results in an overestimate of the contributions to this cell. The
 reverse also applies but there is no guarantee that things cancel out.
 
 Have a look at shade maps of mask and ma0 to see what is happening.
 
 How to fix:
 
 We want an area weighted contribution of all the points on the fine grid.
 
 1) DON'T mask the zeros in the mask variable. We want them for weighting
 correctly.
 2) Use the @AVE regridding to get the correct area weighting.
 
 use mask_15s
 let ma0=mask[gyx=GIN1@AVE]
 list ma0[x=@din,y=@din]/1.e6
 
 should do the trick.
 
 See also
 
 http://ferret.pmel.noaa.gov/Ferret/documentation/users-guide/Grids-Regions/GRIDS#_VPID_198
 
 Cheers,
 Russ
 
 On 25/10/16 23:04, Francois Delclaux wrote:
Hi Ferret users

 I am using a mask at 15" resolution and I calculate the area in km2
 with @DIN operator

 use mask_15s                           !  variable=mask with 0 and 1
 set var/bad=0 mask                       ! 0 is defined as missing value
 list mask[x=@din, y=@din] / 1.E06        ! surface in km2

 the result is correct : 15805. km2


 Now I regrid the mask to 1mn:
 let ma0 = mask[gxy=GIN1]

 and re-compute integral:
 list ma0[x=@din, y=@din] / 1.E06

 the result is: 1.695E-06 , *so the surface is not the same.*

 What happens during the regriding % @DIN operator ?

 Thanks


 *characteristics of variable mask:*
 GRID GAP1
  name axis              # pts   start                end
  LONGITUDE LONGITUDE         1020mr   75.502E(75.502)  79.748E(79.748)
  LATITUDE  LATITUDE           480 r   41.419N              43.415N
  LEVEL2    Z (level)            1 r-  0                    0
  TIME3     T (time)             1 r   0                    0
  normal    E
  normal    F

 stat mask

 MASK
 LONGITUDE: 75.5E(75.5) to 79.8E(79.8)
 LATITUDE: 41.4N to 43.4N
              Z (level): 0
              T (time): 0
              E:  N/A
              F:  N/A
 DATA SET: mask_15s.nc

  Total # of data points: 489600 (1020*480*1*1*1*1)
  # flagged as bad  data: 389847
  Minimum value: 1
  Maximum value: 1
  Mean    value: 1 (unweighted average)
  Standard deviation: 0


 *characteristics of variable ma0*
 GRID GIN1
 name axis # pts start end
 LON X 206 r 75.917 79.333
 LAT Y 76 r 41.75 43
 LEVEL Z 1 r 0 0
 M2 TIME 168 r 16-JAN-2001 05:14 17-DEC-2014 04:14
 normal E
 normal F

 stat ma0
 MASK[GXY=GIN1]
 X: 75.9 to 79.3
 Y: 41.7 to 43
 Z (level): 0
 T (time): 0
 E: N/A
 F: N/A
 DATA SET: mask_15s.nc

 Total # of data points: 15656 (206*76*1*1*1*1)
 # flagged as bad data: 9555
 Minimum value: 1
 Maximum value: 1
 Mean value: 1 (unweighted average)
 Standard deviation: 0


 > Francois DELCLAUX
 > ------------------------------------------------------------
 > UMR HydroSciences - CC 57
 > Université Montpellier
 > 163 rue Auguste Broussonnet
 > 34090  Montpellier      FRANCE
 >
 >http://www.hydrosciences.org/
 > mailto: francois.delclaux@xxxxxxxxxxxxxxxxxxx
 > Tel : (33) (0)4 67 14 90 11      Fax : (33) (0)4 67 14 47 74
 > ------------------------------------------------------------


 --
 François DELCLAUX
 ------------------------------------------------------------
 UMR HydroSciences - CC 57
 Université Montpellier
 163 rue Auguste Broussonnet
 34090 Montpellier FRANCE

 http://www.hydrosciences.org/
 mailto: francois.delclaux@xxxxxxxxxxxxxxxxxxx
 Tel : (33) (0)4 67 14 90 11 Fax : (33) (0)4 67 14 47 74
 ------------------------------------------------------------





[Thread Prev][Thread Next][Index]
Contact Us
Dept of Commerce / NOAA / OAR / PMEL / Ferret

Privacy Policy | Disclaimer | Accessibility Statement