Re: [ferret_users] Problem with external function

[James Caveen <james_caveen@xxxxxxxxxx>]

Hi Ansley,

Thank you for you help.  

I have tried the solution you suggested and indeed it works perfectly
well.  I have asked a colleague of mine to do some further testing with
more voluminous datasets and he should give me some feedback in a few
days.   Meanwhile, I will translate the in-line documentation of my
external function from "la langue de Moliere" to something more
Shakespearean   and as soon as everything is OK, I will forward you the
code.


I have spent a couple of weeks, on and off, trying to make this thing
work and you managed to find a solution within a matter of minutes. So,
once again, thank you for your help.


Regards,

James




Le mercredi 18 avril 2007 à 13:00 -0700, Ansley Manke a écrit :
> 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
> >
> >
> >
> >   
> 
-- 
James Caveen
Analyste-programmeur
Service des technologies de l'information
Université du Québec à Rimouski
310 des Ursulines
Rimouski (Québec) G5L 3A1
Tel: (418) 723-1986 poste 1295
Fax: (418) 724-1842