[Paraview] netCDF reader

[Brockmann Patrick]

>> I have identified those points which causes problems with the actual
>> /ParaView3/VTK/IO/vtkNetCDFCOARDSReader.cxx
>> - if a regular spacing coordinate variable decreases (90 to -90 for
>> example), you always get a vtkRectilinear grid.
>> It should reach, wth the computing choice made, to the creation of
>> vtkImageData. This is the case if you put
>> the axis in an increasing order.
>> You can correct this by placing a fabs set on tolerance variable
>> 108c108
>> < double tolerance = fabs(0.01*this->Spacing);
>> ---
>>     
>>> double tolerance = 0.01*this->Spacing;
>>>       
>
> I'm not sure vtkImageData (and everything that uses it) really supports
> negative spacing.  See, for example, this email thread:
> http://www.vtk.org/pipermail/vtkusers/2009-May/100989.html.
>   
I agree. And that's why I have added the function fabs to make positive the
tolerance variable.

 From ParaView3/VTK/IO/vtkNetCDFCOARDSReader.cxx (line 105...)
    this->Spacing
      = (this->Coordinates->GetValue(dimLen-1) - this->Origin)/(dimLen-1);
    this->HasRegularSpacing = true;     // Then check to see if it is false.
    double tolerance = 0.01*this->Spacing;
    for (size_t i = 1; i < dimLen; i++)
      {
      double expectedValue = this->Origin + i*this->Spacing;
      double actualValue = this->Coordinates->GetValue(i);
      if (   (actualValue < expectedValue-tolerance)
          || (actualValue > expectedValue+tolerance) )
        {
        this->HasRegularSpacing = false;
        fprintf(stderr,"HasRegularSpacing = false");
        break;
        }
      }

As I have written, if a coordinate variable goes from 90 to -90 
(decreasing) with dimLen=7
(-90, -60, -30, 0, 30, 60, 90)
then Spacing variable is (-90 -90)/6 = -30 so negative, as tolerance 
value equals to -0.3
expectedValue for i equals  90 +30 = 120 instead of 60 !
Consequence you get a "HasRegularSpacing = false"  whereas the grid is 
regular.
Put a simple fabs solves the problem
double tolerance = fabs(0.01*this->Spacing);
Is it clearer ?
>   
>> - Test on units should be case independant
>> line 133, units.find(" since ")
>> I have netCDF files where "since" is typed with uppercase (generated by
>> ferret software from the NOAA/PMEL).
>>     
>
> That is easy enough to change.
>   
Do you mean that the user have to change it by himself using ncrename 
(netcdf operators).
I hope not and I consider that it should be treated within the reader. A 
user does not care
about a typology of an attribut.
Having said that, I don't know how to lowercase a vtkStdString.
>> User choices for zorigin and zscale should be taken in consideration .
>> Both set to 1 by default.
>> In addition, a user choice zpositive to tell the direction in which the
>> coordinate values are increasing
>> should also be proposed. With allowed values : up or down (up by default).
>> Perhaps also a zorder choice to treat files where zdimension has to be
>> reordered to get k=1 representing
>> surface.
>>     
>
> I am confused here.  Where are zorigin and zscale defined?  I do not see any
> mention of them in the CF or COARDS specifications.
>   
Indeed CF or COARDS conventions are not concerned by the projection made 
from a lon,lat grid.
But define a height variable as follows
    height = zorigin + fabs(zpositive) * zscale * heightCoords->GetValue(k)
    (with zorigin, zorigin, zscale = 1 by default)
would be better than simply define it as
    height = heightCoords->GetValue(k)
because heightCoords can be negative, have to be scaled, can be equal to 
0 (surface).
>> There is also the problem of connectivity and blank cells (and Z
>> clipping range - in spherical coordinates you
>> don't want to see the globe when missing cells exist).
>>     
>
> Of course I have noticed this.  I have left it this way because it is not
> clear what the right solution is.  There are two issues to closing the grid.
> The first is that the CF/COARDS convention does not really specify whether
> the grid is actually supposed to be closed.  It is perfectly valid to have
> an open sphere (in fact common in the latitudinal dimension).  Of course, a
> simple check of the distance between the two should clear that up.  Or
> perhaps looking at the bounds of the data.
>   
The modulo aspect concerns longitude, not latitude.
In the COARDS convention, one can consider that the bounds of each cell 
is defined at the middle
of 2 center cells.
 From a coordinate variable as follows
lon(imt)= {30, 90, 150, 210, 270, 330}
we can calculate its boundaries
bounds_lon(imt,2) = {0,60, 60,120, 120,180, 180,240, 240,300, 300,360}

There is a modulo longitude to set if first lower cell equals to last 
upper cell (use cosinus).
Here cos(0)=cos(360) so the grid is modulo and the vtkStructuredGrid has 
to be
topologically closed

>> Plus, the fact that CF conventions covers also curvilinear grids where
>> longitudes and latitudes
>> are described with a 2D arrays.
>>     
>
> Uh, isn't this the only way to define curvilinear grids in CF?  It is
> already supported by the netCDF reader.
>   
Not sure of that.
Try from the example file:
http://dods.ipsl.jussieu.fr/prism/gridsCF/sampleCurveGrid4.nc
I get a vtkImageData that represents only the array not its lonlat 
representation.
>> This is why I haved used a vtkPolydata
>> structure (with limitation
>> that the application was not really to visualize 3D field, only 2D
>> fields in a 3D mapping when
>> projection uses spherical coordinates.
>>     
>
> Oh.  I see why you proposed using poly data now.  I still don't think this
> is feasible because of the 3D grid issue.  
However, this what I have done in the python/vtk application
http://www.prism.enes.org/Publications/Reports/Report19.pdf
It uses a vtkPolyData with points defined from the boundaries of each cell.
Either explicitely found from bounds attributs or deducted from middle of
2 center cells.
And indeed a vtkCleanPolyData() filter was passed after creation to get a
clean grid without redondancies.

First goal for me is to get a vtkNetCDFCOARDSReader.cxx working
for grid as defined in the COARDS convention ie a coordinate
variable describes a 1D vector where bounds are at the middle of 2 points.

Patrick