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