[las_users] Troubleshooting Ferret script within LAS

[Tim Cera <timcera@xxxxxxxxxxxxx>]

Trying if visualize a NetCDF file within LAS.  I create the NetCDF file 
from a hydrodynamic model called EFDC.

Here is the dump of the header from my NetCDF file:
----
netcdf efdcout {
dimensions:
        i_rho = 259 ;
        j_rho = 121 ;
        wet_cells = 3013 ;
        nv = 4 ;
        time = UNLIMITED ; // (24 currently)
variables:
        int i_rho(i_rho) ;
        int j_rho(j_rho) ;
        int i_wet_cells(wet_cells) ;
        int j_wet_cells(wet_cells) ;
        float lat_rho(j_rho, i_rho) ;
                lat_rho:long_name = "latitude" ;
                lat_rho:units = "degrees_north" ;
                lat_rho:standard_name = "latitude" ;
                lat_rho:bounds = "lat_rho_bounds" ;
        float lon_rho(j_rho, i_rho) ;
                lon_rho:long_name = "longitude" ;
                lon_rho:units = "degrees_east" ;
                lon_rho:standard_name = "longitude" ;
                lon_rho:bounds = "lon_rho_bounds" ;
        float mask_rho(j_rho, i_rho) ;
                mask_rho:long_name = "mask on RHO-points" ;
                mask_rho:option_0 = "land" ;
                mask_rho:option_1 = "water" ;
                mask_rho:coordinates = "lat_rho lon_rho" ;
                mask_rho:grid_mapping = "transverse_mercator" ;
        float water_elevation(time, j_rho, i_rho) ;
                water_elevation:standard_name = 
"height_above_reference_ellipsoid" ;
                water_elevation:units = "m" ;
                water_elevation:coordinates = "lat_rho lon_rho" ;
                water_elevation:valid_min = -10.f ;
                water_elevation:_FillValue = 1.e+37f ;
        float lat_rho_bounds(nv, j_rho, i_rho) ;
        float lon_rho_bounds(nv, j_rho, i_rho) ;
        int time(time) ;
                time:standard_name = "time" ;
                time:units = "hours since 1995-1-1 0:0:0" ;

// global attributes:
                :Conventions = "CF-1.4" ;
}
----

The cell boundaries of the water cells are 4 sided polygons defined by 
'lat_rho_bounds' and 'lon_rho_bounds'.  The 'lat_rho' and 'lon_rho' 
variables are not used by ferret, but can be used by some NetCDF 
utilities (ncview works great!).

I have a ferret script that visualizes the 'water_elevation' variable 
from the above dataset.  It was hard won and I obtained it by a form of 
genetic algorithm that went something like this, "That didn't work, 
maybe this will..."

Here it is as part of a LAS init script:

----
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! EFDC

! Execute the custom commands to initialize the data. Assumes a dataset and
! variable symbol are defined: data_url and data_var. These are 
unchanged on exit.
! data_var can be the actual variable soil_depth, soil_type, etc.

SET DATA ($data_url)

!<custom

! Define new grid

define axis/x=1:31339:1 vectoraxis
define axis/z=1:4:1 ccaxis
define axis/t=1:24:1 taxis

let llgrid = x[gx=vectoraxis] + z[gz=ccaxis]
let newlat = reshape(lat_rho_bounds, llgrid)
let newlon = reshape(lon_rho_bounds, llgrid)

let flattengrid = x[gx=vectoraxis] + t[gt=taxis]
let newdata = reshape(water_elevation, flattengrid)

polygon/hlimits=-82:-81/vlimits=30.2:30.6 newlon, newlat, newdata[t=1]

!custom>

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
----

Ferret displays what I want, but I have hard coded limits in the polygon 
command and the init script doesn't work anyway in LAS.


Here is a very small part of the LAS error:
----
Server at Revision r3974M
**TMAP ERR: non-existent or not on line
             efdc_polygon_init_script (.jnl)
GO efdc_polygon_init_script 0
Command file, command group, or REPEAT execution aborted
STOP -script mode, ERROR RUNNING SCRIPT statement executed
...
----

Questions:
How do you give the ferret script the horizontal and vertical limits 
within LAS?
Am I just thinking about this all wrong?
Is there an easier way?

Typical curvilinear approach wasn't working within Ferret/LAS because I 
do not have the coordinates for the land cells.  Even if the values are 
masked, seemingly all NetCDF utilities require all cells to have 
coordinates.  I stumbled on the polygon command and not only does it 
visualize what I want, because it is actually using the coordinates of 
each cell corner it honors what the model is using rather than an 
interpolation.

Kindest regards,
Tim Cera