[Thread Prev][Thread Next][Index]

Re: [ferret_users] regriding area larger than grid size

Hi all,
Francois and I looked in to this question a bit, and found the following -

It's true that a SET REGION acts in a highly global way throughout a Ferret session, sometimes doing more than one would expect. Here the setting is interacting with the "subspan modulo" longitude axes that your datasets define to spread the region that the regridding is acting upon beyond the ends of your axis.  This is something we'll look into; the use of subspan modulo axes has resulted in a few other odd interactions.

As far as commenting out the line in your script with define region/x=5:26/y=3:30 tcha ; set reg tcha, you would also need to issue a CANCEL REGION command to actually turn off that setting.  A SET REGION stays in place until you CANCEL REGION, or of course until you exit Ferret and start a new session.


DELCLAUX Francois wrote:
Hi Ansley,

I finaly found the bug !

Your scripts are working fine, so I used them for going back to my script. In fact, at the begining of my script, I put an init command, which defines my interest area, the lake Chad basin. And what happens is that when the following command :
define region/x=5:26/y=3:30 tcha ; set reg tcha
is activated before all the others commands, then I got  wrong results in regriding procedure.
But when I comment this line and start again ferret with my script, it works fine.
So defining a region which does not match exactly the regriding area leads to bad results in regriding.
The second problem is that, even if you suppress the region definition by a comment WHITOUT restarting ferret, the pre-defined region stay active, and so the problem goes on !

Any way, thanks for your help.


François Delclaux

Le 13/10/2009 19:40, Ansley Manke a écrit :
Hello Francois,
Let's look into this, and then we can report the answer back to the Users List when we find out what is happening

I have tried to create an example similar to what you are doing, but I don't see the behavior of the extra points. It's as if Ferret is putting data at the modulo void point  (see "subspan modulo axes" in the documentation) and using that in the regridding step, but this does not happen in the example I created, nor should it happen.

What does the output of ncdump look like for your input files?

  ncdump -c pr_pmip2_0k.nc
  ncdump -c evap_chad.nc   

Here is what I tried, using your variable names so I could be sure to be trying the same commands, but just creating some data.  See what happens if you try these two scripts:
! make_regrid_files.jnl
! Create some variables on the same grids, and write them to files

define axis/units=degrees_east/x=-20:30/npoints=51 xproj
define axis/units=degrees_north/y=-10:30/npoints=41 yproj
define axis/t=1:12:1 tcys

let median = x[gx=xproj] + y[gy=yproj] + 0*t[gt=tcys]
save/clobber/file=pr_pmip2_0k.nc median

define grid/x=xproj/y=yproj/t=tcys glm1
define axis/units=degrees_east/x=7.25:25.25/npoints=37 longitude
define axis/units=degrees_north/y=5.25:23.25/npoints=37 latitude
define axis/z=1:1:1 level
define axis/t=1-jan-1937:09-dec-1938/npoints=708 time

let evap = x[gx=longitude] + y[gy=latitude] + z[gz=level] + 0*t[gt=time]
save/clobber/file=evap_chad.nc evap

Now, in a separate Ferret session, run do_regrid.jnl
! do_regrid.jnl
! Regrid median[d=1] to the grid of evap[d=2]

use pr_pmip2_0k.nc  ! d=1

sh dat
sh grid median

use evap_chad.nc         ! d=2
sh dat
sh grid evap

let p0 = median[d=1, gx=gqw1@lin, gy=gqw1@lin]
save/clobber/file=es.nc p0

DELCLAUX Francois wrote:
Hi Ferret users,

I would like to regrid a GCM output precipitation on a lower  resolution grid defined with another netcdf file. The selected regriding operation is the simplest one, linear interpolation (@lin).

Step one : I read the GCM precip data:
yes? use pr_pmip2_0k.nc  ! d=1

Precip variable name is median and the corresponding grid is GLM1, whose caracteristics are (resol=1°):
yes? sh g glm1
 name       axis              # pts   start                end
 XPROJ     LONGITUDE           51mr   20W(-20)             30E
 YPROJ     LATITUDE            41 r   10S                  30N
 normal    Z
 TCYS      T                   12 r   1                    12

Step 2 : I read the other netcdf file
yes? use evap_chad.nc         ! d=2

Precip variable name is evap and the corresponding grid is GQW1, (resol=0.5°):
yes? sh g gqw1
 name       axis              # pts   start                end
 LONGITUDE LONGITUDE           37mr   7.25E                25.25E
 LATITUDE  LATITUDE            37 r   5.25N                23.25N
 LEVEL     Z (level: "uni       1 r   1                    1
 TIME      TIME               708 r   01-JAN-1937 00:00    09-DEC-1938 00:00

step 3: I regrid median on evap grid  GQW1:
yes? let p0 = median[d=1, gx=gqw1@lin, gy=gqw1@lin]

step 4 : I store p0 in a cdf file and get a look with ncdump :
yes? save/file=es.nc   p0 ; spa ncdump -c es.nc

Result  of dump:
LISTing to file es.nc
netcdf es {
    LONGITUDE0_38 = 39 ;
    LATITUDE = 37 ;
    TCYS = UNLIMITED ; // (12 currently)
    double LONGITUDE0_38(LONGITUDE0_38) ;
        LONGITUDE0_38:units = "degrees_east" ;
        LONGITUDE0_38:modulo = 360. ;
        LONGITUDE0_38:point_spacing = "even" ;
        LONGITUDE0_38:axis = "X" ;
        LATITUDE:units = "degrees_north" ;
        LATITUDE:point_spacing = "even" ;
        LATITUDE:axis = "Y" ;
    double TCYS(TCYS) ;
        TCYS:axis = "T" ;
    float P0(TCYS, LATITUDE, LONGITUDE0_38) ;
        P0:missing_value = -1.e+34f ;
        P0:_FillValue = -1.e+34f ;
        P0:long_name = "MEDIAN[D=1, GX=GQW1@LIN, GY=GQW1@LIN]" ;

// global attributes:
        :history = "FERRET V6.2    8-Oct-09" ;
        :Conventions = "CF-1.0" ;

 LONGITUDE0_38 = -163.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75, 10.25, 10.75,
    11.25, 11.75, 12.25, 12.75, 13.25, 13.75, 14.25, 14.75, 15.25, 15.75,
    16.25, 16.75, 17.25, 17.75, 18.25, 18.75, 19.25, 19.75, 20.25, 20.75,
    21.25, 21.75, 22.25, 22.75, 23.25, 23.75, 24.25, 24.75, 25.25, 196.25 ;

 LATITUDE = 5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
    10.25, 10.75, 11.25, 11.75, 12.25, 12.75, 13.25, 13.75, 14.25, 14.75,
    15.25, 15.75, 16.25, 16.75, 17.25, 17.75, 18.25, 18.75, 19.25, 19.75,
    20.25, 20.75, 21.25, 21.75, 22.25, 22.75, 23.25 ;

 TCYS = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 ;

My problem is that regriding on a 37x37 cell grid (GQW1) creates a 39x37 file , adding for the longitude 2 aditional points on the x-grid longitude : -163.75 and  196.25.

I don't undestand what occurs and how to prevent it.  If someone could help me ...

[Thread Prev][Thread Next][Index]

Contact Us
Dept of Commerce / NOAA / OAR / PMEL / TMAP

Privacy Policy | Disclaimer | Accessibility Statement