Re: [ferret_users] spatial daily mean

[Ansley Manke <ansley.b.manke@xxxxxxxx>]

Hi,
So, you want the average of the data inside the polygonal region for each day.  The function PT_IN_POLY can be used to make a mask. What you need is wht's done in FAQ:
http://ferret.pmel.noaa.gov/Ferret/faq/averaging-integrating-over-irregular-regions/

Here's  an example like the one in the pt_in_poly documentation but for a variable with a time axis

yes? use coads_climatology

yes? ! Define the polygon
yes? let xp = {300,320,335,315,300}
yes? let yp = {20,35,40,20,20}

yes? ! Define the mask
yes? let mask = if pt_in_poly(sst[L=1],xp, yp) ge 0 then 1
yes? let sst_in_poly = mask* sst

The variable "sst_in_poly" has the same data as sst, but is marked as missing everywhere outside the polygon, and it's on the same XYT grid as sst. This means we can take the XY average over the globe or any region that contains the polygon, and the result is the data within the polygon averaged on the grid in XY. 

yes? list sst_in_poly[x=@ave,y=@ave]
             VARIABLE : SST*MASK
             FILENAME : coads_climatology.cdf
             FILEPATH : /home/users/tmap/ferret/linux/fer_dsets/data/
             SUBSET   : 12 points (TIME)
             LONGITUDE: 20E to 20E(380) (XY ave)
             LATITUDE : 90S to 90N (XY ave)
 16-JAN      /  1:  21.80
 15-FEB      /  2:  21.18
 17-MAR      /  3:  21.12
 16-APR      /  4:  21.53
 16-MAY      /  5:  22.55
...

One detail about using this function - you'll notice I defined the mask with a region in L,

   pt_in_poly(sst[L=1],xp, yp)

The polygon definition is in XY and so of course it doesn't depend on the time direction. Why bother to give a time? (By the way any value of L would give exactly the result.)  If I hadn't done that, then the variable "mask" would be 3D - defined on the XY grid and repeated for each value of T. This is no big deal with a short time axis like this, but if there were a long time axis, or a grid with even more dimensions, it's inefficient and would require more memory to load the same mask at each and every time point.

-Ansley

On 5/11/2015 11:20 PM, Nitin Patil wrote:

Dear Ferret users,

In my below script, I have 122 days of rainfall data (l=18050:18171), I have to do the spatial daily mean of rainfall so that I will get 122 mean values in the defined points OR mean of xci,yci grids.

=======================================================================
use rainfall.nc

!! cxi and yci is my region lat lon

let xci = {77,79,79,83,83,84,84,83,83,82,82,80,80,79,79,76,76,73,73,72,72,74,74,75,75,77,77}
let yci = {27,27,25,25,24,24,21,20,20,19,18,18,17,17,16,16,15,15,15,18,21,21,23,23,24,24,24}
let ci = pt_in_poly(rf[d=1],xci, yci)
let ci_data = if ci eq 1 then ci*rf[d=1,l=18050:18171]

list ci_data !! need 122 values listed

=======================================================================

Any suggestion to get it, Kindly let me know.

Thanks,

Nitin