[gdal-dev] WCS GetCoverage with AxisOrder swap

[Even Rouault]

On mercredi 8 novembre 2017 15:50:48 CET Piero Campalani wrote:
> Hi there,
> 
> The gml:GridFunction is *not* an auxiliary definition of the grid geometry,
> but rather a guidance on how to link the grid to the *range set* (i.e. the
> data), so it shall not be compared with the "domain" information (offset
> vectors and so on).
> 
> As Even properly underlines, the order of the coordinates in the offset
> vectors shall refer to the order of the CRS axis definition, that is
> (concerning EPSG:4326) the first coordinate/number is the latitude, the
> second is the longitude (with proper -/+ sign whether the COVERAGE axis
> attached to the vector is concordant/discordant with the CRS axis
> direction).
> 
> The order of the offset-vectors themselves refers to the order of the GRID
> axes (the pic in [1] might be helpful decoding the XML of it all!).
> This order is what can be referred to in the GridFunction: +1 --> first
> GRID axis, +2 --> second GRID axis, etc.
> 
> For instance, what rasdaman seems to tell on the coverage is the the first
> 9000-points COVERAGE axis is the vertical axis, increasing South-ward,
> while the second 18000-points axis follows the East direction.
> Via GridFunction, it then says then that the data you will get in a WCS
> GetCoverage response will span the 2nd (horizontal) axis first.

OK thanks for the explanations. Makes some sense.... 
So it would seem that if you want to return a GeoTIFF in
traditional GIS order (ie faster varying dimension corresponds to longitude), AND you emit
offsetVector in the way Rasdaman does, it seems
that emitting a gml:GridFunction +2 +1 would be required.

But actually, looking closely at the drawing at
http://rasdaman.org/wiki/PetascopeUserGuide#GridaxislabelsandCRSaxislabels , I see
RectifiedGrid.axisLabels="t Long Lat", whereas RectifiedGrid.origin.Point at axisLabels="t Lat Long"
(note the inversion). The second offset vector expressed an increase in longitude (0 0 0.5), while the
third one a decrease in latitude (0 -0.5 0). Which is the GIS friendly order.

So, it would seem that there are 2 more or less equivalent way of achieving the same end result ?

And thus MapServer output at
http://194.66.252.155/cgi-bin/BGS_EMODnet_bathymetry/ows?SERVICE=WCS&REQUEST=DescribeCoverage&version=2.0.1&coverageid=BGS_EMODNET_CentralMed-MCol
could probably be correct, but it would be better for RectifiedGrid.axisLabels to be changed to "long lat" to better
reflect what is done.
Rasdaman output at http://ows.rasdaman.org/rasdaman/ows?SERVICE=WCS&REQUEST=DescribeCoverage&VERSION=2.0.1&COVERAGEID=BlueMarbleCov would be correct
And GeoServer output at https://msp.smartsea.fmi.fi/geoserver/wcs?SERVICE=WCS&REQUEST=DescribeCoverage&VERSION=2.0.1&COVERAGEID=smartsea__eusm2016-EPSG2393
would either need to remove the GridFunction (ala MapServer) or keep it and invert the order in which its offsetVector appear (ala Rasdaman)


And for WCS subsetting, when you write something like SUBSET=AXIS_NAME(min,max),
where AXIS_NAME should come from ? From the RectifiedGrid.axisLabels I guess ?
(to be opposed to CoverageDescription.boundedBy.Envelope at axisLabels)
Then in Ari's test with GeoServer
https://msp.smartsea.fmi.fi/geoserver/wcs?SERVICE=WCS&REQUEST=DescribeCoverage&VERSION=2.0.1&COVERAGEID=smartsea__eusm2016-EPSG2393, in theory "i" and "j" should be the axis requested ?

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20171108/2968671c/attachment.html>