[gdal-dev] Amplitude virtual bands for complex datasets ?

[Even Rouault]

Hi Julien,

> 
> I would like to discuss an idea that I think would be of great interest for
> those using Gdal to access complex dataset, such as SLC SAR products.
> Usually, when one want to visualize such data, looking at the real or
> imaginary part of the signal does not make much sense. People will compute
> amplitude, intensity, or even log-intensity and that is what they will
> look at (phase is to my best knowledge not very useful for visualization,
> but has a lot of value for processing).

Which definition do you give to intensity : amplitude^2 ?

> 
> Computing amplitude and phase from the complex value could be left to the
> end-user software, but : - Most software do not provide such features
> (think of Qgis for instance, or MapServer), - If overviews are generated
> from the real/imaginary image values with interpolation (nearest neighbor
> would be fine), then the generated overviews will be wrong, because you
> cannot interpolate complex data this way.

AFAICS in the code, GDAL currently implements for complex data types : nearest 
neighbour, average on both real and imaginary terms, and AVERAGE_MAGPHASE 
which is average but with renormalization of the amplitude so that the  
AVERAGE_MAGPHASE'ed amplitude is the average of the source amplitudes.

> 
> I think that Gdal could greatly ease the pain by exposing to the user
> virtual subdatasets corresponding to the amplitude, phase, intensity and
> log-intensity, for all complex datasets. In the case of multi-band complex
> datasets, those virtual subdatasets should also be multi-band, so that
> polarimetric data could be accessed directly as a multi-band amplitude
> image for instance.

One of my fear for exposing such computed subdatasets is that there could be a 
confusion between products that are made of "natural" subdatasets (several 
images of different size/resolution in the same container), and those computed 
subdatasets. For example if you use gdal_translate -sds, you probably only 
want "natural" subdatasets to be converted.
Could also get tricky to implement for drivers that can support "natural" 
subdatasets and complex data. Sentinel1 SLC products for example.

In your idea would the generation of those virtual subdatasets would/could be 
done in GDAL core, or should that be enabled per-driver ? Because of the 
example I gave with Sentinel1, I'm not sure a completely generic solution can 
be found. Perhaps something mix with a generic part & per-driver adaptations 
to "merge" natural and derived subdatasets.

Another solution would be to use VRT derived bands (see "Using Derived Bands" 
in http://www.gdal.org/gdal_vrttut.html) and provide a few hard-coded pixel 
functions to compute amplitude, phase, etc...
and have API facilities like
GDALRasterBandH GDALCreateDerivedBand( GDALRasterBandH hSrCband, const char* 
pszPixelFunction [, GDALDataType eTargetDT ?] ) that would create the in-
memory VRT. But I'm aware that doesn't address your idea of using unmodified 
QGIS, MapServer, etc...

> 
> From what I know from Gdal capabilities this is not much work. I could try
> to do it myself, but I would like to get your feedback first (and some
> hints on where to start would be useful to me).

There would be likely changes to do in:
- GDALOpen(), to detect a syntax like 
"DERIVED_SUBDATASET:Amplitude:original_datasetname" and create the appropriate 
dataset. Hum, or perhaps better, instead of hacking GDALOpen(), having a 
driver that would recognize this syntax and instanciate the proper derived 
dataset.
- GDALDataset::GetMetadata( pszDomain ) when pszDomain equals "SUBDATASETS"
-> with the issue I mentionned for drivers that already implement subdatasets 
and will not call the base method.


Hum, or perhaps a better solution would be to have a new metadata domain 
"DERIVED_SUBDATASETS". Would solve the potential confusion between "natural" 
and derived subdatasets, and would probably require little/no changes in 
drivers. Would not confuse gdal_translate -sds. Could be used by mapserver 
unmodified. But would require QGIS querying this new metadata domain.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com