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.
Ansley
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.
Sincerely
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
GRID 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
GRID 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 {
dimensions:
LONGITUDE0_38 = 39 ;
LATITUDE = 37 ;
TCYS = UNLIMITED ; // (12 currently)
variables:
double LONGITUDE0_38(LONGITUDE0_38) ;
LONGITUDE0_38:units = "degrees_east" ;
LONGITUDE0_38:modulo = 360. ;
LONGITUDE0_38:point_spacing = "even" ;
LONGITUDE0_38:axis = "X" ;
double LATITUDE(LATITUDE) ;
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" ;
data:
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 ...
|