[gdal-dev] WCS GetCoverage with AxisOrder swap

Piero Campalani piero.campa at gmail.com
Wed Nov 8 06:50:48 PST 2017


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.

hth!
-Piero

[1]
http://rasdaman.org/wiki/PetascopeUserGuide#GridaxislabelsandCRSaxislabels
[2]
http://rasdaman.org/wiki/PetascopeSubsets#Geometricinterpretationofacoverage

On 8 November 2017 at 14:42, <gdal-dev-request at lists.osgeo.org> wrote:

> Message: 3
> Date: Wed, 8 Nov 2017 13:48:49 +0200
> From: Ari Jolma <ari.jolma at gmail.com>
> To: "'gdal-dev at lists.osgeo.org' (gdal-dev at lists.osgeo.org)"
>         <gdal-dev at lists.osgeo.org>
> Subject: [gdal-dev] WCS GetCoverage with AxisOrder swap
> Message-ID: <fff6a965-8ac8-ebbc-c244-f81451fc3117 at gmail.com>
> Content-Type: text/plain; charset=utf-8; format=flowed
>
> I've got to the point in the new WCS driver development, where I'm
> trying to issue GetCoverage request to a GeoServer, which is serving
> data in EPSG that has axis order swapped.
>
> This is related to the earlier discussion
>
> https://lists.osgeo.org/pipermail/gdal-dev/2017-April/046366.html
>
> The data serving is at
>
> https://msp.smartsea.fmi.fi/geoserver/wcs?SERVICE=WCS&
> REQUEST=DescribeCoverage&VERSION=2.0.1&COVERAGEID=
> smartsea__eusm2016-EPSG2393
>
> The boundedBy tells that the axisLabels are "Y X" and the lower corner
> is "6543350.381089335 3061592.6391462097" and upper corner is
> "7307056.899896299 3432141.361396798". This makes sense, in that EPSG X
> (easting) runs up from 3000000 and Y (northing) runs up from 6500000.
>
> What I don't understand is why I need to make the GetCoverage request
> with subset=Y(3061592,3061632) and subset=X(7307016,7307056). If I
> define them as I would expect
>
> https://msp.smartsea.fmi.fi/geoserver/wcs?SERVICE=WCS&
> REQUEST=GetCoverage&VERSION=2.0.1&COVERAGEID=smartsea__
> eusm2016-EPSG2393&FORMAT=image%2Ftiff&SUBSET=X%283061592.63914621,3061632.
> 65520693%29&SUBSET=Y%287307016.88383558,7307056.8998963%29
>
> I get "Empty intersection after subsetting".
>
> A MapServer that serves data with EPSG that has axis order swapped is 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
>
> It works using the above logic. However, it does not have GridFunction,
> which GeoServer has (it is "+2 +1", i.e., the data is east first, then
> north).
>
> and Rasdaman server at
> ​​
>
>
> http://ows.rasdaman.org/rasdaa32c78aman/ows?SERVICE=WCS&
> REQUEST=DescribeCoverage&VERSION=2.0.1&COVERAGEID=BlueMarbleCov
> <http://ows.rasdaman.org/rasdaman/ows?SERVICE=WCS&REQUEST=DescribeCoverage&VERSION=2.0.1&COVERAGEID=BlueMarbleCov>
>
> Works using the above logic. However, it *has* GridFunction, which is
> the same as GeoServer has. However, it has the order of offsetVectors
> swapped - which I think is correct (as is also written here
> https://lists.osgeo.org/pipermail/gdal-dev/2017-April/046421.html)
>
> So, I can make Rasdaman work for CRS with swapped order by default (also
> the returned geotiff is ok), and for MapServer and GeoServer I need to
> have at least two hack options (NoOffsetSwap and NoGridEnvelopeSwap,
> maybe there could be only one). I have not yet checked if the rasters
> returned from those two are ok.
>
> Ari
>
>
>
>
>
>
> ------------------------------
>
> Message: 4
> Date: Wed, 8 Nov 2017 14:50:54 +0200
> From: Ari Jolma <ari.jolma at gmail.com>
> To: "'gdal-dev at lists.osgeo.org' (gdal-dev at lists.osgeo.org)"
>         <gdal-dev at lists.osgeo.org>
> Subject: Re: [gdal-dev] WCS GetCoverage with AxisOrder swap
> Message-ID: <bd7723d6-9fd7-fa2c-4770-7c68d4c357d8 at gmail.com>
> Content-Type: text/plain; charset=utf-8; format=flowed
>
> Ari Jolma kirjoitti 08.11.2017 klo 13:48:
> >  for MapServer and GeoServer I need to have at least two hack options
> > (NoOffsetSwap and NoGridEnvelopeSwap, maybe there could be only one).
> > I have not yet checked if the rasters returned from those two are ok.
>
> The rasters that those servers return are ok. I just need the one
> additional hack option for GeoServer to swap the axis names for the
> GetCoverage request.
>
> I also don't need to use the GridFunction information.
>
> BTW, it seems that GDAL driver options are usually all upcase. Is that
> just a convention, and could it perhaps be violated? The reason for
> options like "PreferredFormat=xxx" for the WCS driver (instead of
> "PREFERREDFORMAT=xxx") is that 'PreferredFormat' is the name of the
> element in the WCS_GDAL XML. My approach in the new WCS driver is to use
> options to set values in the WCS_GDAL XML since that XML file is in the
> cache and meant to be maintained by the driver.
>
> Ari
>
>
>
> ------------------------------
>
> Message: 5
> Date: Wed, 08 Nov 2017 14:42:51 +0100
> From: Even Rouault <even.rouault at spatialys.com>
> To: gdal-dev at lists.osgeo.org
> Subject: Re: [gdal-dev] WCS GetCoverage with AxisOrder swap
> Message-ID: <2058494.1KC9oNcdcJ at even-i700>
> Content-Type: text/plain; charset="iso-8859-1"
>
> On mercredi 8 novembre 2017 14:50:54 CET Ari Jolma wrote:
> > Ari Jolma kirjoitti 08.11.2017 klo 13:48:
> > >  for MapServer and GeoServer I need to have at least two hack options
> > > (NoOffsetSwap and NoGridEnvelopeSwap, maybe there could be only one).
> > > I have not yet checked if the rasters returned from those two are ok.
> >
> > The rasters that those servers return are ok. I just need the one
> > additional hack option for GeoServer to swap the axis names for the
> > GetCoverage request.
> >
> > I also don't need to use the GridFunction information.
>
> One thing that is confusing, and rings some bells to me as I saw that with
> GMLJP2 (the
> GDAL_JP2K_ALT_OFFSETVECTOR_ORDER config option in gdaljp2metadata.cpp) is
> the order of the offetVector element
>
> >From your example, we have:
>
> MapServer:
> <gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/0/4326">0
> 0.004167</gml:offsetVector>
> <gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/0/4326">-0.004167
> 0</gml:offsetVector>
>
> GeoServer:
> <gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/0/2393">0.0
> 20.00803035910304</gml:offsetVector>
> <gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/0/2393
> ">-20.00803035910304 0.0</gml:offsetVector>
>
> Rasdaman
> <offsetVector srsName="http://ows.rasdaman.org/def/crs/EPSG/0/4326">-0.02
> 0</offsetVector>
> <offsetVector srsName="http://ows.rasdaman.org/def/crs/EPSG/0/4326">0
> 0.02</offsetVector>
>
> So in all 3 cases, the order of the value *inside* a offetVector element
> is lat/northing long/easting as those CRS as lat-long / northing-easting. OK
>
> BUT Rasdaman does an extra swap of the offsetVector elements themselves.
>
> What is not clear (among many things) is when GridFunction applies. If it
> applies before linking the raster/grid space to the georeferenced space,
> that is the first offsetVector applies to the fastest varying dimension of
> the final grid,
> then I'd say MapServer & GeoServer do the right thing (GeoServer probably
> closer to correctenss, due to te
> presence of the GridFunction). Otherwise, if its applied at the very end
> of the process, Rasdaman is probably correct.
>
> Practically I'd say if you have
> <offsetVector srsName="...">A B</offsetVector>
> <offsetVector srsName="...">C D</offsetVector>
> with SRS with lat-long or northing-easting order, with high confidence
> regarding how "correct" an implementation is:
>  * if A != 0 and D > 0 and C=B=0, D the x resolution and A is the y
> resolution .
>  * if A=D=0 and C != 0 and B > 0, then B is the x resolution and C the y
> resolution.
> In other situations, good luck...
>
> --
> 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/9a786fbe/attachment.html>
>
> ------------------------------
>
> Subject: Digest Footer
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> ------------------------------
>
> End of gdal-dev Digest, Vol 162, Issue 22
> *****************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20171108/ba1c5b3d/attachment-0001.html>


More information about the gdal-dev mailing list