WCPS should automatically detect 2D transposed data and apply a transpose operation

[Vlad Merticariu]

Overview:

rendering 2D results should follow Long - Lat axis ordering, as this is what the clients and the computer screen follow for maps / images. The current behavior of rasdaman is to follow the order of axis of the array on disk. rasdaman has a transpose operation that can be used to deal with such cases. Detecting Lat - Long cases can be done automatically, and a transpose operation can be injected in the WCPS query.

Details / test cases:

There are cases where the Long - Lat ordering is respected, even if the CRS of a coverage says it is Lat - Long, e.g:

​http://ows.rasdaman.org/rasdaman/ows?&SERVICE=WCS&VERSION=2.0.1&REQUEST=DescribeCoverage&COVERAGEID=RadianceColor

​http://ows.rasdaman.org/rasdaman/ows?&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=RadianceColor&SUBSET=ansi(%222002-07-01T00:00:00.000Z%22)&FORMAT=image/png

This happens because the real sdom of the RadianceColor object is [0:3599,0:1799,0:149], so on disk the data is Long, Lat.

However, in the case below, both the crs and the data on disk are Lat, Long, thus the 2D result is Lat Long too.

​http://earthserver.pml.ac.uk/rasdaman/ows?&SERVICE=WCS&VERSION=2.0.1&REQUEST=DescribeCoverage&COVERAGEID=OCCCI_V3_monthly_rrs_443_bias

​http://earthserver.pml.ac.uk/rasdaman/ows?&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=OCCCI_V3_monthly_rrs_443_bias&SUBSET=ansi(%221997-09-04T00:00:00.000Z%22)&FORMAT=image/png

Detecting the cases where the returned image could work in the following way:

• check if the result is 2D
• if yes, check the real order of axis in the returned data. This can be done by sorting the vector of Geo axes according to their rasdamanOrder property.
• if the resulting vector is Lat Long (or better Y X in axisType), add "{ \"transpose\": [0,1] }" to the WCPS query.

[Peter Baumann]

I'm not in favour of doing this automatically and invisibly. Rather, users could be provided with a WCPS function transpose() which operates on the axis labels.

[Vlad Merticariu]

We already offer that, unfortunately this doesn't work with WCS.
This is strictly a usability decision., affecting only the visualization of the data.

For example, requesting the above mentioned data through WCS for display on a WorldWind globe doesn't work (​​http://earthserver.pml.ac.uk/rasdaman/ows?&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=OCCCI_V3_monthly_rrs_443_bias&SUBSET=ansi(%221997-09-04T00:00:00.000Z%22)&FORMAT=image/png).

[Peter Baumann]

got it. So what about defining a transposition request parameter as a facade for the corresponding WCPS transpoe() ?

[Dimitar Misev]

I think the proper way would be to transpose at ingest time. But I guess that won't generally work well for insitu as it will be a waste to transpose if we don't export a lat/lon result.

[Vlad Merticariu]

I don't think we should transpose at import, as we want to keep the data as provided by the user.

This is just a heuristic determining that the result is intended to be displayed in a X / Y order.

I am in favor of doing this automatically because that is the behavior of similar tools:

• gdal automatically transposes the data if translated from a netcdf variable to a GeoTiff for example
• geoserver automatically transposes netcdf data when displayed as png

[Dimitar Misev]

Ok I think it makes sense in petascope as we surely know what the axis order of Lat/Lon should be. Does this heuristic cover all cases of CRS, e.g. axes E/W as well?

[Vlad Merticariu]

In theory it should, as we follow axis directions (x and Y), but we're at the mercy of the CRS correctly defining these. For the existing ones that's covered, but we can't control future ones.

[Peter Baumann]

indeed, the CRS definition is what we must rely on. Simply swapping axes would invalidate the coverage as axes would not match with the CRS. What we still could do is a different sequence rule, changing "+1 +2" to "-2 +1" or so.
But I concur with Dimitar: what is the general rule? In particular, for 3D+?

[Vlad Merticariu]

In the gmlcov model, the connection between the CRS axis sequence and the data axis sequence is given by the offset vectors. In particular, a coverage having offset vectors:

• 0 1
• 1 0

indicates that the first data axis corresponds to the second axis of the CRS, and the second data axis corresponds to the first axis of the CRS.
Transposing the data would mean changing the offset vectors to:

• 1 0
• 0 1

with the first data axis corresponding now to the first CRS axis.

[Vlad Merticariu]

The general rule on how we store the data is: same ordering as the data used to create the coverage. Where new dimensions are added to the data (e.g. creating a cube out of 2D files), the extra dimensions can be added left or right (e.g. time first or time last). Serving the data with any other axis ordering (e.g. always creating a 1 to 1 match between CRS ordering and data ordering) could be done by transposing, but it would be too costly to do for all queries. Doing it at insert time works, but only for imported data, not for insitu.

Given that we want to transpose as little as possible, we know when the user wants it to happen when the format is for direct displayed of 2D maps, and the data axis ordering is not X Y, so we could do it just in these cases.

We can identify that the 2D result of a WCS / WCPS query has data axes Y X by one of:

• maintaining a list of axisLabels that are generally X (Long, Longitude, E, Y etc.) and one for Y (Lat, Latitude, N, etc.)
• check the <gml:axisDirection> element in the CRS body, and hope that all Lat axes point "north", which I know is not really reliable

Currently we implement the first option, but that is fragile to custom axis labels. We could improve this by allowing the user to indicate that an axis is of type X.

The alternative to this, mentioned earlier, would be to always serve data with the axis order indicated by the CRS. Would this be better from a standards' perspective? Right now, the CRS alone is not enough to know how to take some data and put it on a map. But the axis ordering is not the only element that prevents you from doing that, you also need the pixel size (which is also encoded in the offset vector).