[gdal-dev] Multidimensional raster support in GDAL

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] Multidimensional raster support in GDAL

Ari Jolma-2
I'd like to get back to this topic a bit since I learned something from
the standards that I believe was not discussed in the earlier discussion

https://lists.osgeo.org/pipermail/gdal-dev/2017-October/047464.html

In WCS 1 the range of the data - remember the main conceptual division
to the domain (the X,Y, etc dimensions) and the range (the data if you
will) - consists of fields, which may consist of bands. That is, a field
may have an array value (an array of bands). In WCS 2, this separation
is, firstly, outsourced to SWE since the datatype of range is
SWE::DataRecord, and, secondly, left as further work in the range
subsetting extension. The current range subsetting extension makes no
difference between fields and bands, they are the same and there is no
structure.

The GDAL WCS driver has had some experimental support for time as field
and different time steps as bands. I'm not sure if there are many people
using it since I believe the functionality is currently broken when used
together with MapServer. The brokenness seems to be so serious that when
making a gdal_translate call where the source is a WCS that has a field
with more than one band (and GDAL recognizes it as such) and the call
does not involve scaling, the call fails. The reason seems to be that
gdal_translate accesses the data band by band.

I guess my point is here is that this is another difference between the
GDAL data model and the general coverage data model. The simple approach
to the data model mismatch is to allow GDAL datasets to be made from
coverage sources with queries. Also one thing I realized I need for the
new WCS driver is something which distinguishes different GDAL datasets
made from a single coverage.

Best regards,

Ari


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Multidimensional raster support in GDAL

Even Rouault-2

On mardi 21 novembre 2017 12:23:31 CET Ari Jolma wrote:

> I'd like to get back to this topic a bit since I learned something from

> the standards that I believe was not discussed in the earlier discussion

>

> https://lists.osgeo.org/pipermail/gdal-dev/2017-October/047464.html

>

> In WCS 1 the range of the data - remember the main conceptual division

> to the domain (the X,Y, etc dimensions) and the range (the data if you

> will) - consists of fields, which may consist of bands. That is, a field

> may have an array value (an array of bands). In WCS 2, this separation

> is, firstly, outsourced to SWE since the datatype of range is

> SWE::DataRecord, and, secondly, left as further work in the range

> subsetting extension. The current range subsetting extension makes no

> difference between fields and bands, they are the same and there is no

> structure.

 

Looking at the WCS 1.1 spec, the RangeSubset KVP syntax shows indeed a lot of complexity.

 

RangeSubset = FieldSubset *( “;” FieldSubset )

FieldSubset = FieldId [ “:” Interpolation ] [ “[” AxisSubsets “]” ] AxisSubsets = AxisSubset *( “,” AxisSubset )

AxisSubset = AxisId “[” Keys “]”

Keys = Key *( “,” Key )

 

Not sure if implemetations in practice go beyond having a single field made of a single axis (at least this is what MapServer supports).

 

>

> The GDAL WCS driver has had some experimental support for time as field

 

I'd say more than the time support is a kind of ad-hoc dimensional parameter.

 

> and different time steps as bands.

 

My reading of the code and doc of the current driver is that time steps are treated as subdatasets, and not bands.

See https://trac.osgeo.org/gdal/ticket/3449

 

> I'm not sure if there are many people

> using it since I believe the functionality is currently broken when used

> together with MapServer. The brokenness seems to be so serious that when

> making a gdal_translate call where the source is a WCS that has a field

> with more than one band (and GDAL recognizes it as such) and the call

> does not involve scaling, the call fails.

 

Is it really an issue specific of time or something more general about multiband support ?

 

I've just spent some time looking at mapserver code & doc, and there's some confusion regarding multiband support . In particular the WCS 2.0 code has a mode where it uses the WCS 1.x mapfile metadata item "wcs_rangeset_axes" in a way not compatible with the WCS 1.x implementation, see

https://github.com/mapserver/mapserver/pull/5523

 

However for WCS 1.0, I managed to make band subsetting work with the following conf in the mapfile (no time component)

 

"wcs_bandcount" "9"

"wcs_rangeset_axes" "bands"

"wcs_rangeset_name" "Landsat 5 TM Bands"

"wcs_rangeset_label" "Bands"

"wcs_rangeset_description" "Bands for Landsat 5 TM"

 

DescribeCoverage reports:

 

<rangeSet>

<RangeSet>

<description>Bands for Landsat 5 TM</description>

<name>Landsat 5 TM Bands</name>

<label>Bands</label>

<axisDescription>

<AxisDescription>

<name>bands</name>

<label>Bands/Channels/Samples</label>

<values>

<singleValue>1</singleValue>

<singleValue>2</singleValue>

<singleValue>3</singleValue>

<singleValue>4</singleValue>

<singleValue>5</singleValue>

<singleValue>6</singleValue>

<singleValue>7</singleValue>

<singleValue>8</singleValue>

<singleValue>9</singleValue>

</values>

</AxisDescription>

</axisDescription>

<nullValues>

<singleValue>0</singleValue>

</nullValues>

</RangeSet>

</rangeSet>

 

And the following GetCoverage request succeeds and returns 3 bands

SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&COVERAGE=multi&FORMAT=GEOTIFF_8&BBOX=15,48,16,49&bands=9,5,1&CRS=EPSG:4326&WIDTH=5&HEIGHT=5

 

> The reason seems to be that

> gdal_translate accesses the data band by band.

 

gdal_translate (actually GDALDatasetCopyWholeRaster()) uses the INTERLEAVE=PIXEL/BAND metadata item of the IMAGE_STRUCTURE domain as a hint whether to read band by band or all bands together

 

So if you add the following, it should request all bands :

poDS->SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE")

 

>

> I guess my point is here is that this is another difference between the

> GDAL data model and the general coverage data model. The simple approach

> to the data model mismatch is to allow GDAL datasets to be made from

> coverage sources with queries. Also one thing I realized I need for the

> new WCS driver is something which distinguishes different GDAL datasets

> made from a single coverage.

 

I'm not sure what you have in mind,

 

>

> Best regards,

>

> Ari

>

>

> _______________________________________________

> gdal-dev mailing list

> [hidden email]

> https://lists.osgeo.org/mailman/listinfo/gdal-dev

 

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Multidimensional raster support in GDAL

Ari Jolma-2
Even Rouault kirjoitti 21.11.2017 klo 16:06:

 

> and different time steps as bands.

 

My reading of the code and doc of the current driver is that time steps are treated as subdatasets, and not bands.


Yes. My mistake.

 

> I'm not sure if there are many people

> using it since I believe the functionality is currently broken when used

> together with MapServer. The brokenness seems to be so serious that when

> making a gdal_translate call where the source is a WCS that has a field

> with more than one band (and GDAL recognizes it as such) and the call

> does not involve scaling, the call fails.

 

Is it really an issue specific of time or something more general about multiband support ?

 

I've just spent some time looking at mapserver code & doc, and there's some confusion regarding multiband support . In particular the WCS 2.0 code has a mode where it uses the WCS 1.x mapfile metadata item "wcs_rangeset_axes" in a way not compatible with the WCS 1.x implementation, see

https://github.com/mapserver/mapserver/pull/5523

 

However for WCS 1.0, I managed to make band subsetting work with the following conf in the mapfile (no time component)

 

"wcs_bandcount" "9"

"wcs_rangeset_axes" "bands"

"wcs_rangeset_name" "Landsat 5 TM Bands"

"wcs_rangeset_label" "Bands"

"wcs_rangeset_description" "Bands for Landsat 5 TM"

 

DescribeCoverage reports:

 

<rangeSet>

<RangeSet>

<description>Bands for Landsat 5 TM</description>

<name>Landsat 5 TM Bands</name>

<label>Bands</label>

<axisDescription>

<AxisDescription>

<name>bands</name>

<label>Bands/Channels/Samples</label>

<values>

<singleValue>1</singleValue>

<singleValue>2</singleValue>

<singleValue>3</singleValue>

<singleValue>4</singleValue>

<singleValue>5</singleValue>

<singleValue>6</singleValue>

<singleValue>7</singleValue>

<singleValue>8</singleValue>

<singleValue>9</singleValue>

</values>

</AxisDescription>

</axisDescription>

<nullValues>

<singleValue>0</singleValue>

</nullValues>

</RangeSet>

</rangeSet>

 

And the following GetCoverage request succeeds and returns 3 bands

SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&COVERAGE=multi&FORMAT=GEOTIFF_8&BBOX=15,48,16,49&bands=9,5,1&CRS=EPSG:4326&WIDTH=5&HEIGHT=5


So it's a misconfiguration at 194.66.252.155

http://194.66.252.155/cgi-bin/BGS_EMODnet_bathymetry/ows?SERVICE=WCS&REQUEST=GetCoverage&VERSION=1.0.0&COVERAGE=BGS_EMODNET_CentralMed-MCol&FORMAT=GTiff&BBOX=9.83125,44.05624829,14.09792008,46.18958333&WIDTH=1024&HEIGHT=512&CRS=EPSG:4326&bands=1

returns 3 bands although the DescribeCoverage is similar.

 

> The reason seems to be that

> gdal_translate accesses the data band by band.

 

gdal_translate (actually GDALDatasetCopyWholeRaster()) uses the INTERLEAVE=PIXEL/BAND metadata item of the IMAGE_STRUCTURE domain as a hint whether to read band by band or all bands together

 

So if you add the following, it should request all bands :

poDS->SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE")


ok, with that I can remove the need for the user option. And probably it is anyway usually better,

 

>

> I guess my point is here is that this is another difference between the

> GDAL data model and the general coverage data model. The simple approach

> to the data model mismatch is to allow GDAL datasets to be made from

> coverage sources with queries. Also one thing I realized I need for the

> new WCS driver is something which distinguishes different GDAL datasets

> made from a single coverage.

 

I'm not sure what you have in mind,


Somehow the user must be able to write

dataset_1 = GDAL.Open("WCS:URL?coverage=coverage_a")
dataset_2 = GDAL.Open("WCS:URL?coverage=coverage_a")

So that for example dataset_1 is band c from coverage_a and dataset_2 is band d from the same coverage (this example is stupid but I think this is a valid issue with more complex coverages). I guess that needs to be in the dataset name somehow since that is the key to to service file.

Ari


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Multidimensional raster support in GDAL

Even Rouault-2

> > And the following GetCoverage request succeeds and returns 3 bands

> >

> > SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&COVERAGE=multi&FORMAT=GEOTIF

> > F_8&BBOX=15,48,16,49&bands=9,5,1&CRS=EPSG:4326&WIDTH=5&HEIGHT=5

> So it's a misconfiguration at 194.66.252.155

 

Yes, probably an issue with the definition of their GTiff output format

See

https://github.com/mapserver/mapserver/issues/3919#issuecomment-4944131

 

> Somehow the user must be able to write

>

> dataset_1 = GDAL.Open("WCS:URL?coverage=coverage_a")

> dataset_2 = GDAL.Open("WCS:URL?coverage=coverage_a")

>

> So that for example dataset_1 is band c from coverage_a and dataset_2 is

> band d from the same coverage (this example is stupid but I think this

> is a valid issue with more complex coverages). I guess that needs to be

> in the dataset name somehow since that is the key to to service file.

 

You probably intended to add band specification in the above dataset names ?

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Multidimensional raster support in GDAL

Ari Jolma-2
Even Rouault kirjoitti 21.11.2017 klo 18:51:

 

> Somehow the user must be able to write

>

> dataset_1 = GDAL.Open("WCS:URL?coverage=coverage_a")

> dataset_2 = GDAL.Open("WCS:URL?coverage=coverage_a")

>

> So that for example dataset_1 is band c from coverage_a and dataset_2 is

> band d from the same coverage (this example is stupid but I think this

> is a valid issue with more complex coverages). I guess that needs to be

> in the dataset name somehow since that is the key to to service file.

 

You probably intended to add band specification in the above dataset names ?


I intentionally left them out, since up to now I've used the URL with version and coverage name as the cache key but it's clear that all parameters that can affect the contents of the dataset must be in the name (URL, cache key) and not in options.

Ari


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Multidimensional raster support in GDAL

Ari Jolma-2
In reply to this post by Even Rouault-2
Even Rouault kirjoitti 21.11.2017 klo 16:06:

 

gdal_translate (actually GDALDatasetCopyWholeRaster()) uses the INTERLEAVE=PIXEL/BAND metadata item of the IMAGE_STRUCTURE domain as a hint whether to read band by band or all bands together

 

So if you add the following, it should request all bands :

poDS->SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE")


I don't see that happening in the case of WCS driver:

PIXEL_INTERLEAVE=TRUE is my option to set the metadata item but it could be set with option INTERLEAVE=PIXEL, which GDALDatasetCopyWholeRaster() examines.

8<-------------------
gdal_translate -oo PIXEL_INTERLEAVE=TRUE -oo OriginAtBoundary -oo CACHE=wcs_cache -srcwin 0 0 2 2 "WCS:http://194.66.252.155/cgi-bin/BGS_EMODnet_bathymetry/ows?version=1.0.0&coverage=BGS_EMODNET_CentralMed-MCol" x
GNM: GNMRegisterAllInternal
GNM: RegisterGNMFile
GNM: RegisterGNMdatabase
Setting INTERLEAVE to PIXEL
GDAL: GDALOpen(WCS:http://194.66.252.155/cgi-bin/BGS_EMODnet_bathymetry/ows?version=1.0.0&coverage=BGS_EMODNET_CentralMed-MCol, this=0x1256210) succeeds as WCS.
Input file size is 2977, 3883
0GTiff: Reopen with strip chop enabled
GDAL: GDAL_CACHEMAX = 585 MB
GDAL: GDALDatasetCopyWholeRaster(): 2*2 swaths, bInterleave=1
URL=http://194.66.252.155/cgi-bin/BGS_EMODnet_bathymetry/ows?SERVICE=WCS&REQUEST=GetCoverage&VERSION=1.0.0&COVERAGE=BGS_EMODNET_CentralMed-MCol&FORMAT=GTiff&BBOX=9.83125,44.05624829,14.09792008,46.18958333&WIDTH=1024&HEIGHT=512&CRS=EPSG:4326&bands=1
HTTP: Fetch(http://194.66.252.155/cgi-bin/BGS_EMODnet_bathymetry/ows?SERVICE=WCS&REQUEST=GetCoverage&VERSION=1.0.0&COVERAGE=BGS_EMODNET_CentralMed-MCol&FORMAT=GTiff&BBOX=9.83125,44.05624829,14.09792008,46.18958333&WIDTH=1024&HEIGHT=512&CRS=EPSG:4326&bands=1)
HTTP: libcurl/7.47.0 OpenSSL/1.0.2g zlib/1.2.8 libidn/1.32 librtmp/2.3
WCS: GDALOpenResult() on content-type: image/tiff
GDAL: GDALOpen(/vsimem/wcs/0x1256210/wcsresult.dat, this=0x12a2610) succeeds as GTiff.
ERROR 1: Returned tile does not match expected band configuration.
Got 3 bands instead of one although the coverage has band range type.

GDAL: GDALClose(/vsimem/wcs/0x1256210/wcsresult.dat, this=0x12a2610)
ERROR 1: wcs_cache/znQnc.xml, band 1: IReadBlock failed at X offset 0, Y offset 0: Returned tile does not match expected band configuration.
Got 3 bands instead of one although the coverage has band range type.

GDAL: GDALClose(wcs_cache/znQnc.xml, this=0x1256210)
GDAL: In GDALDestroy - unloading GDAL shared library.
8<-------------------



Ari


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Multidimensional raster support in GDAL

Even Rouault-2

On mercredi 22 novembre 2017 14:38:55 CET Ari Jolma wrote:

> Even Rouault kirjoitti 21.11.2017 klo 16:06:

> > gdal_translate (actually GDALDatasetCopyWholeRaster()) uses the

> > INTERLEAVE=PIXEL/BAND metadata item of the IMAGE_STRUCTURE domain as a

> > hint whether to read band by band or all bands together

> >

> > So if you add the following, it should request all bands :

> >

> > poDS->SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE")

>

> I don't see that happening in the case of WCS driver:

>

> PIXEL_INTERLEAVE=TRUE is my option to set the metadata item but it could

> be set with option INTERLEAVE=PIXEL, which GDALDatasetCopyWholeRaster()

> examines.

>

 

Probably related to the logic of WCSDataset::TestUseBlockIO() which given your area of interest size of 2x2 will force bocked IO instead of direct IO

 

If PIXEL_INTERLEAVE=TRUE is set, you might want to force the implementation of WCSRasterBand::IReadBlock() to read all the bands.

 

Note: the above is with the current trunk logic. Didn't check if you modified this.

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Multidimensional raster support in GDAL

Ari Jolma-2
Even Rouault kirjoitti 22.11.2017 klo 15:11:

On mercredi 22 novembre 2017 14:38:55 CET Ari Jolma wrote:

> Even Rouault kirjoitti 21.11.2017 klo 16:06:

> > gdal_translate (actually GDALDatasetCopyWholeRaster()) uses the

> > INTERLEAVE=PIXEL/BAND metadata item of the IMAGE_STRUCTURE domain as a

> > hint whether to read band by band or all bands together

> >

> > So if you add the following, it should request all bands :

> >

> > poDS->SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE")

>

> I don't see that happening in the case of WCS driver:

>

> PIXEL_INTERLEAVE=TRUE is my option to set the metadata item but it could

> be set with option INTERLEAVE=PIXEL, which GDALDatasetCopyWholeRaster()

> examines.

>

 

Probably related to the logic of WCSDataset::TestUseBlockIO() which given your area of interest size of 2x2 will force bocked IO instead of direct IO

 

If PIXEL_INTERLEAVE=TRUE is set, you might want to force the implementation of WCSRasterBand::IReadBlock() to read all the bands.


I guess the situation is same for WCSDataset::IRasterIO which may be called BandCount and BandMap.

Works fine.

Ari


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev