[gdal-dev] WCS GetCoverage with AxisOrder swap

Thu Nov 9 01:53:54 PST 2017

Even,

On 8 November 2017 at 16:31, Even Rouault <even.rouault at spatialys.com>
wrote:

> 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.
>

Yes, as Jukka also pointed out in a later reply.  !

If you /match/ the description of the geometry of the grid (origin, offset
vectors order, etc) with the way you list the grid _data_ then you don't
need the GridFunction (that is, you can assume the default +1 +2 ... +N
linear rule).

In the other case, you need to instruct the client on how to interpret the
data.

(I believe that if the GridFunction would not be optional, this all would
be a little less confusing)

> But actually, looking closely at the drawing at
>
> http://rasdaman.org/wiki/PetascopeUserGuid
> <http://rasdaman.org/wiki/PetascopeUserGuide#GridaxislabelsandCRSaxislabels>
> 
> <http://rasdaman.org/wiki/PetascopeUserGuide#GridaxislabelsandCRSaxislabels>
> e#GridaxislabelsandCRSaxislabels
> <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 thethird 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 ?
>

GRID axes names (in spite of the CRS axes names which are formally
described in its definition) are arbitrary: in the example, the rasdaman
implementation chooses to name grid axes by the CRS axes they span.

RectifiedGrid.axisLabels (GRID axes) are "t Long Lat", but could just be
"0 1 2", or "a b c". Naming them after the CRS they span is just more
convenient.

RectifiedGrid.origin.Point at axisLabels are the CRS axes names instead,
hence latitude goes first.
The GML there is OK: indeed you see the 2nd offset-vector, the direction of
the 2nd "Long"-labelled GRID axis) goes '0 0 0.5', hence in the easting
direction.

>
> And thus MapServer output at
>
> http://194.66.252.1
> <http://194.66.252.155/cgi-bin/BGS_EMODnet_bathymetry/ows?SERVICE=WCS&REQUEST=DescribeCoverage&version=2.0.1&coverageid=BGS_EMODNET_CentralMed-MCol>
>
> 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.
>
>
Exactly, that output is not OK as it is now.

> 
> 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)
>
>
>

If GeoServer is not going to span the North direction first in the
GetCoverage responses, then that description is wrong.
Being EPSG:2393 a North-first CRS, the geometry of the grid already defines
the North-ward as first (+1), and East-ward as second (+2), so if that is
the rule which GeoServer will use to list the data, the GridFunction there
is wrong.

> 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 ?
>

NO, subsetting/subspacing is meant to work aligned with the coordinate
system (!)
You choose your (sub)space of interest, and then the server will filter out
all the grid points that get included in that space, independently of its,
say, geometric profile (indeed in extreme settings, what you get might not
be a grid! see below)

Example: imagine an oblique 100x100 2D regular grid in a 3D "X Y Z" CS
space. The grid has 2 "i" and "j" axes, with origin in "0 0 0" and offset
vectors which might be "1 1 0.5" for i and "-0.5 -0.5 1" for j.
The user's might want to fetch all the grid points in the small cube
delimited by X(0:3) Y(0:3) Z(7:10). in the 3D space, the user has defined
this little cube which will act as spatial filter for the oblique grid.

Subsetting via GRID indexes would be probably achieved if you request a
GetCoverage by using the "Index?D" CRS which is mean

Anyway.. it was just for keeping it a little less abstract.
WCS subsettings need the CRS axis names.

(In this nasty example, if you'd want to subset the oblique grid by taking
its 10x10 subgrid starting from the origin, then a) the server would need
to implement the WCS CRS extension, b) you would need to set the
"subsettingCrs" parameter to, say,
http://ows.rasdaman.org/def/crs/OGC/0/Index2D in the GetCoverage request,
and subset to i(0:9) and j(0:9))
[1] http://rasdaman.org/wiki/IndexCrss

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 ?
>
>
No: having explained that here above, you see subsets in Ari's case should
be "X" and "Y", as defined in EPSG:2393.

-
Piero !
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20171109/0fb4d88f/attachment.html>