NetCDF positive attribute in vertical Z axis variable seems to be ignored

[SMichalet]

The netCDF attribute "positive" for vertical Z axis seems to be ignored. I take a netCDF file whith "positive" attribute equal to "down" for a vertical Z Axis (depth in my netcdf file) :

float depth(depth) ;
        depth:long_name = "depth" ;
        depth:units = "m" ;
        depth:positive = "down" ;

with values :

depth = 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 35, 40,
45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 100, 125, 150, 175, 200, 250,
300, 350, 400, 500, 600, 700, 800, 900, 1000, 1250, 1500, 2000, 2500,
3000, 4000, 5000 ;

And I get positive values in ncWMS getCapabilities :

<Dimension name="elevation" units="m" default="0.0">
0.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 100.0, 125.0, 150.0, 175.0, 200.0, 250.0, 300.0, 350.0, 400.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0, 1250.0, 1500.0, 2000.0, 2500.0, 3000.0, 4000.0, 5000.0
</Dimension>

In the first version of ncWMS, with the same file, I get correct negative values :

<Dimension name="elevation" units="m" default="-0.0">
-0.0,-2.0,-4.0,-6.0,-8.0,-10.0,-12.0,-14.0,-16.0,-18.0,-20.0,-22.0,-24.0,-26.0,-28.0,-30.0,-35.0,-40.0,-45.0,-50.0,-55.0,-60.0,-65.0,-70.0,-75.0,-80.0,-85.0,-90.0,-100.0,-125.0,-150.0,-175.0,-200.0,-250.0,-300.0,-350.0,-400.0,-500.0,-600.0,-700.0,-800.0,-900.0,-1000.0,-1250.0,-1500.0,-2000.0,-2500.0,-3000.0,-4000.0,-5000.0
</Dimension>

I'm using ncWMS 2.2.3. Am I missing something ?

[guygriffiths]

That looks like a bug. I'll try and find some time to fix it, but I'm currently very busy on other projects at the moment.

[guygriffiths]

In the libraries which underpin ncWMS2, we keep metadata about the vertical CRS (including the fact that positive values increase downwards). The problem is that this metadata is not making it through to the WMS capabilities document - there is no way of distinguishing whether a WMS layer's vertical coordinates represent depth or height, which client applications may want for layer comparison.

However, I don't think there's anything intrinsically wrong with the capabilities document stating that the values of depth on the vertical axis are positive (since they are...). So I think that the underlying issue is that we're not specifying the units correctly - reading the WMS spec, units should indicate the CRS of the vertical data, and the unitSymbol attribute should reflect the units of measurement within that vertical CRS.

If the vertical axis values represent height in metres, we can use units="CRS:88" as recommended in the WMS spec. If they represent depth in metres, then either units="EPSG:5715" or units="EPSG:5831" would be a suitable choice (depth from MSL and instantaneous sea-level respectively, see http://epsg.io/5715 and http://epsg.io/5831). If they are not measured in metres, then technically we should be using a CRS whose units match. However, there aren't any named ones in the EPSG database, so I think that we should just bend the WMS spec and use the same CRS for units + the appropriate units for unitsSymbol - it's going to be very rare that data specifies depth/height in units of length other than metres, and we're already not following the spec entirely correctly.

The problem comes when we have units of pressure, where there are no named vertical CRSs defined. The best option in that case may be to use an ad hoc approach, perhaps defining the units attribute as barometric and the unitsSymbol as whatever the pressure is measured in.

@jonblower - What do you think about this?

[jonblower]

there is no way of distinguishing whether a WMS layer's vertical coordinates represent depth or height

I think in ncWMS1 the vertical axis was always "elevation" but if the underlying data was "depth" we just used negative values in the Capabilities doc. Could this work in EDAL/ncWMS2? Using the CRS codes is possible but I don't know how many clients will recognise them.

For pressure I can't remember what we did, I suspect we just set the units to "Pa" or something like that.

[guygriffiths]

The problem with just negating the values in the Capabilities document is that GetMap will need to negate requested elevation values as well, which would mean that we would have to do that throughout EDAL. So for an axis with a vertical CRS stating that positive is down and containing values of 0,5,10, our general-purpose methods for finding whether a value is contained in the axis would have to report true for -5 and -10 and false for 5 and 10. We could do that and document the behaviour, but it's very unintuitive and is bound to cause headaches.

The only way to feasibly do it is at the data reading level - we'd need to detect the positive="down" attribute, make the vertical CRS state that positive is up and then negate the values on axis creation. At which point we should remove the isPositiveUpwards() metadata from the VerticalCrs class (since it will only be false for pressure-based systems, for which we already have isPressure()). But I'm not really sure that that wouldn't cause issues with dimensionless quantities such as sigma levels.

Either way, we should probably use CRS codes for the units attribute in vertical <dimension> tags in the Capabilities document, since we're not currently following the spec. I'm not sure that it's important that clients actually recognise CRS codes - I think that the important thing is whether the vertical units of 2 different layers are the same, and hence can be directly compared, but I guess it depends on the use case.

@SMichalet - What issues is this bug causing you? Knowing that would help a lot in deciding what the best solution is.

[jonblower]

Yes, I think if we simply negate the values we'd only do this "at the last minute" in the creation of the Capabilities doc, or somewhere in the WMS-specific code. In EDAL the metadata should be as close to the data files as possible (i.e. usually with positive coordinate values and appropriate return values for isPositiveUpwards().

IIRC, in ncWMS1, we did something like "if isPositiveUpwards() == false and axisType != pressure then use an "elevation" axis but negate the coordinate values". (I might be wrong though.)

One thing to watch is that giving a CRS code might imply a level of specificity that we can't warrant from the data (EPSG 5715 and 5831, and "depth below ellipsoid" are all different CRSs and we can't always tell the difference from the source files). I don't know what the best "default" behaviour might be.

[guygriffiths]

So no-one from the ncwms-users list seems to have an opinion, and after talking to Martin Desruisseaux it seems that the issue of vertical CRSes is not settled amongst either the OGC or in the CF-conventions.

As we discussed offline, I think the best solution is that we use the units attribute to store a CRS-like string, specifically either ncwms:height, ncwms:depth, or ncwms:pressure and then store the actual units in the unitSymbol attribute. This would move us closer to WMS compliance, and whilst the "CRSes" we would be using would just be arbitrary strings, they do classify the vertical CRS into the 3 categories which can currently be distinguished by the CF conventions.