Re: [ferret_users] Problem with external function

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

Hi James,
Your external function sounds like a useful one; when you're happy with 
it might I ask you to post it back here, or send it to us at 
contact_ferret@xxxxxxxxxxxxx and we can make it available to users?

When you "set data filename" in Ferret, Ferret reads the header 
information from the file and saves it, and also reads and stores 
coordinate information. If an axis has the attribute point_spacing= 
"even" then Ferret doesn't even read the whole list of coordinates; it 
just gets the first and last coordinates, computes the delta-coordinate 
value, and stops there, as that's all that is needed to define a 
regularly-spaced axis. If there is no point_spacing attribute, or if 
it's marked as uneven, then Ferret reads all the coordinates.  It checks 
to see if they are in fact even and if so stores only the first, last, 
and delta.

Ferret does this checking no matter what, since there's no guarantee 
that an axis with the same name from a new file will be the same as one 
in an already-open file, and so in your case the code sees that the new 
axis with a large number of irregular coordinates exceeds the memory 
available for coordinate storage.

Here is a way to force Ferret to use the IWET axis from the first file 
for other files:

1) Open the first file with the correct, irregular IWET axis.
2) Change the attribute for IWET in the other files to point_spacing= 
"even".
  When opened, the grid for variables in these files will have a different
  name for that axis, like IWET1 or IWET2 and Ferret will not read the 
coordinates.
3) Regrid the variable in file 2 to IWET using the @ASN regridding 
transformation.

yes? use file_1
yes? use file_2
yes? let/D=2 temp_2 = temp[gt=iwet@asn]

I haven't experimented with this, but it should work fine.

To keep the name of the variable intact,  you can use SET VAR/NAME= to 
rename the variable after opening the file, then use the original name 
in the regridding definition.

  yes? use file_2
  yes? set var/name=temp_reg temp   ! Give the file variable a new name
  yes? let/D=2 temp = temp_reg[gt=iwet@asn]

Ansley

James Caveen wrote:
> Dear ferreters,
>
> I have been trying to solve a problem with an external function for
> quite some time now and I finally decide to turn to the Ferret
> Community.  I realize that  this message is somewhat long but here I go
> anyway.
>
>  
> I have written a Ferret  external function to reconstruct a 3D grid that
> is stored using compression by gathering as described in Section 8.2 of
> the CF Convention :
> (http://www.cgd.ucar.edu/cms/eaton/cf-metadata/CF-current.html#gath)
>
> and everything works fine when I use one NetCDF file.  However, if I try
> to open a second file, Ferret's limit for irregular axis is exceeded;
> my data is stored in a 295000 element vector with uneven spacing. 
>
> Following is a ncdum -h of one of my files:
>
>
> netcdf COMP_TEST.ROM_TEMP_copie1 {
> dimensions:
>         XGULF = 236 ;
>         YGULF = 150 ;
>         ZGULF = 73 ;
>         ZGULFedges = 74 ;
>         TGULF = UNLIMITED ; // (4 currently)
>         IWET = 295000 ;
> variables:
>         float XGULF(XGULF) ;
>                 XGULF:point_spacing = "even" ;
>         float YGULF(YGULF) ;
>                 YGULF:point_spacing = "even" ;
>         float ZGULF(ZGULF) ;
>                 ZGULF:positive = "down" ;
>                 ZGULF:point_spacing = "uneven" ;
>                 ZGULF:edges = "ZGULFedges" ;
>         float ZGULFedges(ZGULFedges) ;
>         float TGULF(TGULF) ;
>                 TGULF:units = "HOURS since 1901-01-15 00:00:00" ;
>                 TGULF:time_origin = "15-JAN-1901" ;
>         int IWET(IWET) ;
>                 IWET:compress = "XGULF YGULF ZGULF" ;
>                 IWET:point_spacing = "uneven" ;
>         float TEMP(TGULF, IWET) ;
>                 TEMP:fill_value = 9.96921e+36f ;
>                 TEMP:missing_value = 9.96921e+36f ;
>                 TEMP:long_name = "TEMP" ;
>                 TEMP:units = "celcius" ;
>
> and here is an example script where I use my nc_rar() external function:
>
> set memory/size=500
> use COMP_TEST.ROM_TEMP.nc
> let a = x[gx=xgulf]*0 + y[gy=ygulf]*0 + z[gz=zgulf]*0 + t[gt=tgulf]*0
> let b = nc_rar(a,temp])
> shade b[l=1,k=1]
>
>         
>
> So here I go with my Question 1:
> All my files use the same IWET axis.  Is there a way that I can tell
> Ferret that the IWET axis of file two is identical to the IWET axis of
> file one and thus force Ferret not to allocate extra space for the
> second IWET axis?
>
>
> Being somewhat naive, I attempted to fool Ferret by pretending the IWET
> axis is evenly spaced as follows:
>
>  
>         int IWET(IWET) ;
>                 IWET:compress = "XGULF YGULF ZGULF" ;
>                 IWET:point_spacing = "even" ;
>
> Now I can open two,three,... NetCDF files without encountering the
> memory limit problem but my external function doesn't work anymore (all
> indices returned by ef_get_coordinates() are wrong).  
>
> So, here is my second question:
> Could someone briefly explain to me what pre-treatment Ferret does to
> axis variables before passing them on to an external function?
>
>
> Finally, being somewhat stubborn, I decided to convert my IWET axis into
> a regular variable (WET) and to define an evenly spaced axis XWET which
> is common to all my variables:   
>
> netcdf COMP_TEST.ROM_TEMP_copie1 {
> dimensions:
>         XGULF = 236 ;
>         YGULF = 150 ;
>         ZGULF = 73 ;
>         ZGULFedges = 74 ;
>         TGULF = UNLIMITED ; // (4 currently)
>         XWET = 295000 ;
> variables:
>         float XGULF(XGULF) ;
>                 XGULF:point_spacing = "even" ;
>         float YGULF(YGULF) ;
>                 YGULF:point_spacing = "even" ;
>         float ZGULF(ZGULF) ;
>                 ZGULF:positive = "down" ;
>                 ZGULF:point_spacing = "uneven" ;
>                 ZGULF:edges = "ZGULFedges" ;
>         float ZGULFedges(ZGULFedges) ;
>         float TGULF(TGULF) ;
>                 TGULF:units = "HOURS since 1901-01-15 00:00:00" ;
>                 TGULF:time_origin = "15-JAN-1901" ;
>         int XWET(XWET) ;
>                 IWET:compress = "XGULF YGULF ZGULF" ;
>                 IWET:point_spacing = "even" ;
>         float TEMP(TGULF, XWET) ;
>                 TEMP:fill_value = 9.96921e+36f ;
>                 TEMP:missing_value = 9.96921e+36f ;
>                 TEMP:long_name = "TEMP" ;
>                 TEMP:units = "celcius" ;
>         float WET(XWET) ;
>                 WET:fill_value = 9.96921e+36f ;
>                 WET:missing_value = 9.96921e+36f ;
>                 WET:long_name = "Indices points mouilles" ;
>
>
> I then modified my external function so as to supply to it the WET
> variable which contains the wet-cell indices:
>
> set memory/size=500
> use COMP_TEST.ROM_TEMP.nc
> let a = x[gx=xgulf]*0 + y[gy=ygulf]*0 + z[gz=zgulf]*0 + t[gt=tgulf]*0
> let b = nc_rar(a,temp,wet)
> shade b[l=1,k=1]
>
>
> The problem is that now inside my nc_rar() function all variables (temp
> and wet) have zero values everywhere.  Everything else inside my
> function appears to be right: arg_lo_ss, arg_hi_ss, etc. so I obviously
> have an addressing problem or something similar. 
>
> So here is my final question:
>
> Has anybody ever encountered this behaviour with external functions and
> managed to solve it?
>
>
> I realize that Ferret does not support the CF Convention for NetCDF
> files but the disk space savings resulting from compression by gathering
> are very interesting.  So, I intend to investigate this avenue further
> for our own needs.  Any tip or advice to help me solve this problem
> would be greatly appreciated.
>
>
>
>
> James
>
>
>
>   

--

~>   ~>   ~>   ~>   ~>   ~>   ~>   ~>   ~>   ~>   ~> 

Ansley Manke, NOAA/PMEL   ansley.b.manke@xxxxxxxx
7600 Sand Point Way NE, Seattle WA
Phone 206-526-6246,  FAX 206-526-6744