[Thread Prev][Thread Next][Index]

Re: set region subtlety (and point correlation script)

Hello everyone,

I want to add a further clarification of the question Rob brought up.
The heart of the question is in the application of a region when using a
grid-changing function such as RESHAPE.  See this FAQ for more on
SET REGION and grid-changing functions:

A note right off -- Rob asked if it takes more memory to define a region
with a qualifier on the variable, as in:

   let sst_np = sst[GX=xax_np,y=35N:45N@AVE]

The answer is that it does not take more memory, and often it is the best
or only way to be clear about the region you want to apply to a variable.

If LIST commands are inserted into Rob's script at the different stages
of the calculation, it gives correct results up until the RESHAPE operation.
Now, RESHAPE is a grid-changing function.  Any limits made with
SET REGION apply to the result grid of the expression involving
RESHAPE, not to the input arguments.  The arguments need their
own region specifications, when a region is needed.

When sst_np is defined with qualifying region definition,

        let sst_np = sst[GX=xax_np,y=35N:45N@AVE]

then the RESHAPE operation at the end of the script has as its first
argument a variable that contains its own region information, Y=35:45.
This Y region is applied to the variable before it is used in RESHAPE.

But in the version of the script with

        set region /y=35N:45N;  let sst_np = sst[GX=xax_np,y=@AVE]

the first argument to RESHAPE does not contain its own Y limits AND
IT WILL NOT inherit the limits given in SET REGION.   Again, the
SET REGION limits apply to the *result* grid of the expression, which
is the grid of "dummy" in this script, and dummy has no Y limits.  Since
no Y limits are given within the sst_np definition, and all the subsequent
variables that depend on sst_np, the evaluation of RESHAPE will
default to the full span of the sst Y axis.  (One should be able to see
this difference using SET MODE DIAGNOSTIC output.)  To correctly
impose the Y limits you could insert "Y=35:45" to modify the first argument:

   reshape(covar[y=35N:45N], dummy[k=4])

OR, best of all, just use the definition where sst_np has the region qualifier

    sst_np = sst[GX=xax_np,y=35N:45N@AVE]

Ansley Manke

Rob Scott wrote:

Ferret Users:

I've encountered a problem using set region that makes me uneasy.
Should these two lines,

        set region /y=35N:45N
        let sst_np = sst[GX=xax_np,y=@AVE]

not give the same numerical results as these:

        let sst_np = sst[GX=xax_np,y=35N:45N@AVE]

?? I avoided the latter because I thought it was harder on memory, no?

Well they give different results in the following script

Rob Scott

(hopefully someone may find it useful for doing point correlations)

! Description: get the point correlation of 1D SST field

set data coads_climatology
define axis/x=140E:130W:1/units=longitude/edges xax_np

! the following gets the wrong answer

set region /y=35N:45N
let sst_np = sst[GX=xax_np,y=@AVE]

! but this would make it work:
!            **     let sst_np = sst[GX=xax_np,y=35N:45N@AVE]   **

! keep q fixed:
! move p in zonal direction, find acvf = <p q>
let q = sst_np[x=145W]
! setup NetCDF file to hold the answer
let dummy = x[GX=xax_np]
let/title="SST Autocovariance" acvf = (1/0) * x[GX=xax_np]
save/clobber/file=sst_acvf.cdf acvf
! for all p points, get acvf
! save to NetCDF (use RESHAPE to get it on the right grid)
! save to ASCII (for comparison)
repeat/k=1:5 (let p = sst_np[i=`k`]; go variance; \
LIST/FORMAT=(1X,E14.7)/FILE="sst_acvf.dat"/CLOBBER/NOHEAD covar ; \
LET acvf = reshape(covar,dummy[i=`k`]) ; \
SAVE/APPEND/FILE=sst_acvf.cdf acvf)

Postal Address:

Program in Atmospheric and Oceanic Sciences
P.O. Box CN710, Sayre Hall
Princeton, NJ 08544-0710

Tel: 609-452-6519                    o__    ____ 
Fax: 609-987-5063                    _,>/'_  -----
E-mail: rscott@princeton.edu        (_) \(_) ------

Ansley Manke  Pacific Marine Environmental Laboratory  Seattle WA  (206)526-6246

[Thread Prev][Thread Next][Index]

Dept of Commerce / NOAA / OAR / PMEL / TMAP

Contact Us | Privacy Policy | Disclaimer | Accessibility Statement