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

classic Classic list List threaded Threaded
25 messages Options
12
Reply | Threaded
Open this post in threaded view
|

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

Julien Michel-2
Dear Gdal developers,

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).

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.

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.

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).

Regards,

Julien

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

Re: Amplitude virtual bands for complex datasets ?

Even Rouault-2
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
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Amplitude virtual bands for complex datasets ?

Antonio Valentino
Hi Even, hi Julien,

Il 28/06/2016 12:39, Even Rouault ha scritto:

> 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.

[CUT]


> 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...
>

Sometime ago I proposed to include in gdal a small set of pixel
function.  My use case at the time was mainly relater to SAR SLC data
visualization.

There is an open ticket #3367 [1] an also a git repo with pixel
functions implemented as gdal plugin [2]


[1] https://trac.osgeo.org/gdal/ticket/3367
[2] https://github.com/avalentino/gdal-pixfun-plugin

I hope this can help is some way.
It would be very nice to have an improved support for SAR complex data.

>> 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
>

ciao

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

Re: Amplitude virtual bands for complex datasets ?

Julien Michel-2
Hi Even, Hi Antonio,

The VRT with custom PixelFunction was actually the first thing I thought of. However, if I understand well those function need to be registered to gdal before opening the VRT, and this must be done in client code, so it will not allow support in Qgis or MapServer for instance. Of course, if Gdal could ship a set of standard pixel functions, the VRT solution would be a solid first step. Maybe this is simpler than derived datasets ? I still think the latter would be nicer, even if I understand that Gdal must avoid mixing natural and derived subdatatsets.

Thanks a lot for your answers,

Julien



-----Message d'origine-----
De : gdal-dev [mailto:[hidden email]] De la part de Antonio Valentino
Envoyé : mardi 28 juin 2016 22:50
À : [hidden email]
Objet : Re: [gdal-dev] Amplitude virtual bands for complex datasets ?

Hi Even, hi Julien,

Il 28/06/2016 12:39, Even Rouault ha scritto:

> 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.

[CUT]


> 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...
>

Sometime ago I proposed to include in gdal a small set of pixel function.  My use case at the time was mainly relater to SAR SLC data visualization.

There is an open ticket #3367 [1] an also a git repo with pixel functions implemented as gdal plugin [2]


[1] https://trac.osgeo.org/gdal/ticket/3367
[2] https://github.com/avalentino/gdal-pixfun-plugin

I hope this can help is some way.
It would be very nice to have an improved support for SAR complex data.

>> 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
>

ciao

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

Re: Amplitude virtual bands for complex datasets ?

Bugbuster
In reply to this post by Antonio Valentino
Dear Antonio, dear Even

I have used the PIXFUN plugin a while ago under Linux, and as far as I remember, it worked pretty well  .
I have downloaded the last version a few days ago because I need to run some tests under Windows. I had to add the export instruction for the DLL compilation (CPL_DLL macro defined in gdal_port.h). The plugin compiles correctly, it is registered at startup.
However, it seems not to work correctly. I tested the real(), imag() and mod() functions. Each time, the number of source rasters (nSources) is zero (I added a few log to find this out)
Do you have any idea what could cause this weird behavior ?

Thanks
Philippe
Reply | Threaded
Open this post in threaded view
|

Re: Amplitude virtual bands for complex datasets ?

Antonio Valentino
Hi Bugbuster,

Il 30/06/2016 17:16, Bugbuster ha scritto:

> Dear Antonio, dear Even
>
> I have used the PIXFUN plugin a while ago under Linux, and as far as I
> remember, it *worked pretty well*  .
> I have downloaded the last version a few days ago because I need to run some
> tests under *Windows*. I had to add the export instruction for the DLL
> compilation (CPL_DLL macro defined in gdal_port.h). The plugin compiles
> correctly, it is registered at startup.
> However, it seems not to work correctly. I tested the real(), imag() and
> mod() functions. Each time, the number of source rasters (/nSources/) is
> zero (I added a few log to find this out)
> Do you have any idea what could cause this weird behavior ?

the package on githib also contains a test suite
Can you please add a test that triggers the issue on windows?
I will try to look into it.

Also the patch with the CPL_DLL would be appreciated


cheers

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

Re: Amplitude virtual bands for complex datasets ?

Bugbuster
Hi Antonio,

I did try the test suite you provided. I tried to convert all VRT images to GTiff with the simple command :
gdal_translate -of GTiff pixfun_xx.vrt pixfun_xx.vrt.tif

None of them gave any results.
In each function, I added a debug message when testing the number of source rasters :
if (nSources !=1) {
   CPLDebug("PIXFUN", "ImagePixelFunc() : nSources = %d\n", nSources);
   return CE_Failure);
}
I always get :
nSources = 0

If you have any ideas, I would be very grateful to you :-)
Thanks
Philippe
Reply | Threaded
Open this post in threaded view
|

Re: Amplitude virtual bands for complex datasets ?

Antonio Valentino
Hi Bugbuster,

Il 01/07/2016 08:52, Bugbuster ha scritto:
> Hi Antonio,
>
> I did try the test suite you provided. I tried to convert all VRT images to
> GTiff with the simple command :
>
>
> None of them gave any results.

please ensure to have gdal_PIXFUN.so in the GDAL_DRIVER_PATH:

env GDAL_DRIVER_PATH=../../.. gdal_translate -of GTiff pixfun_cmul_c.vrt
pixfun_cmul_c.vrt.tif

> In each function, I added a debug message when testing the number of source
> rasters :
>
> I always get :
>
>
> If you have any ideas, I would be very grateful to you :-)
> Thanks
> Philippe


regards

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

Re: Amplitude virtual bands for complex datasets ?

Bugbuster
Hi Antonio,

as I was mentioning it in my 1st message,I am running under Windows.
I can see that the PIXFUN is indeed recognized thanks to the log messages I have added in all functions.
I have noticed that the pointer papoSources is always null !

If you have a clue, I would buy it :-)
Sincerly
Philippe
Reply | Threaded
Open this post in threaded view
|

Re: Amplitude virtual bands for complex datasets ?

Even Rouault-2
Le lundi 04 juillet 2016 09:30:44, Bugbuster a écrit :
> Hi Antonio,
>
> as I was mentioning it in my 1st message,I am running under *Windows*.
> I can see that the PIXFUN is indeed recognized thanks to the log messages I
> have added in all functions.
> I have noticed that the pointer *papoSources* is always *null *!
>
> If you have a clue, I would buy it :-)

Hi,

I've add a look at the code and noticed on unsafe practice, but I'm not sure
it explains what you observe. In pixelfunctions.c,
GDALRegisterDefaultPixelFunc() is declared with CPL_STDCALL convention, but in
pixfunplugin.c with "extern CPLErr GDALRegisterDefaultPixelFunc()"; (and thus
C convention I guess).

Even

> Sincerly
> Philippe
>
>
>
>
> --
> View this message in context:
> http://osgeo-org.1560.x6.nabble.com/gdal-dev-Amplitude-virtual-bands-for-c
> omplex-datasets-tp5273713p5274571.html Sent from the GDAL - Dev mailing
> list archive at Nabble.com.
> _______________________________________________
> gdal-dev mailing list
> [hidden email]
> http://lists.osgeo.org/mailman/listinfo/gdal-dev

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Amplitude virtual bands for complex datasets ?

Julien Michel-2
In reply to this post by Even Rouault-2
Le 28/06/2016 à 12:39, Even Rouault a écrit :

> 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.
To check if I understood well :

I will create a driver that will recognize the
"DERIVED_SUBDATASET:Amplitude:original_datasetname" syntax. This driver
needs to now if "original_datasetname" is of complex type (to report it
can open it), and will also expose the derived band somehow (can I
benefit from the CreateDerivedBand API, and the pixel functions from
Antonio ?). Also, "original_datasetname" can be a subdataset itself.

There should also be a mechanism that will report the existing
DERIVED_SUBDATASETS upon query (this I did not get yet).

I am not sure if I can figure out how to do this, but I will try.

Regards,

Julien

--
Julien MICHEL
CNES - DCT/SI/AP

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

Re: Amplitude virtual bands for complex datasets ?

Even Rouault-2
Hi,
>
> To check if I understood well :
>
> I will create a driver that will recognize the
> "DERIVED_SUBDATASET:Amplitude:original_datasetname" syntax. This driver
> needs to now if "original_datasetname" is of complex type (to report it
> can open it)

That or we could also possibly accept non-complex types (assume imaginary part
= 0). Antonio's functions handle complex & non-complex AFAICS

> , and will also expose the derived band somehow (can I
> benefit from the CreateDerivedBand API,

If you write it, yes... (it doesn't exist yet if I wasn't clear)

If you follow the VRT pixel function road, that's a matter of building a on-
the-fly VRT dataset.

http://www.gdal.org/gdal_vrttut.html#gdal_vrttut_creation could possibly be
used. But AFAICS setting the pixel function name isn't supported. Could be
done by implementing CPLErr VRTSourcedRasterBand::SetMetadataItem( const char
*pszName, const char *pszValue, const char *pszDomain ), similarly to what
VRTSourcedRasterBand::SetMetadataItem(...) handles for vrt sources.

> and the pixel functions from Antonio ?).

Would make sense to internalize this code (uses X/MIT license)

> Also, "original_datasetname" can be a subdataset itself.

Yes, that shouldn't matter to the driver. It is just a string that it provides
to GDALOpen() (some care must be taken when parsing
DERIVED_SUBDATASET:algname:original_datasetname, since original_datasetname
can have columns too)

>
> There should also be a mechanism that will report the existing
> DERIVED_SUBDATASETS upon query (this I did not get yet).

Implement in GDALDataset, the following method (from GDALMajorObject)

virtual char **GetMetadata( const char * pszDomain = "" );

char **GDALDataset::GetMetadata( const char * pszDomain = "" )
{
        if( pszDomain != NULL && EQUAL(pszDomain, "DERIVED_SUBDATASET") )
        {
                // some stuff. Be careful that ownership of returned char** belongs
                // to the dataset object
        }
        else
                return GDALMajorObject::GetMetadata(pszDomain);
}

As well as :

virtual char **GetMetadataDomainList(); to report DERIVED_SUBDATASET when it
makes sense (contrary to GetMetadata(), ownership of the returned char**
belongs to the caller)

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Amplitude virtual bands for complex datasets ?

Julien Michel-2
Hi,

Started the work here :

https://github.com/OSGeo/gdal/compare/trunk...jmichel-otb:enhance-complex-datasets

For now, recognizes the following syntax :

gdal_translate -srcwin 1000 1000 1000 1000
DERIVED_SUBDATASET:COMPLEX_AMPLITUDE:s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff
~/tmp/test.tif

I need to see about exposing those derived datasets properly.

Antonio, I copied one of your pixel function for fast prototyping, I
hope you do not mind. Of course upon completion of this work we should
merge properly all your pixel functions with proper
credits/copyright/licence. Thing is I did not know where to put them
within gdal.

Regards,

Julien

Le 04/07/2016 à 11:06, Even Rouault a écrit :

> Hi,
>> To check if I understood well :
>>
>> I will create a driver that will recognize the
>> "DERIVED_SUBDATASET:Amplitude:original_datasetname" syntax. This driver
>> needs to now if "original_datasetname" is of complex type (to report it
>> can open it)
> That or we could also possibly accept non-complex types (assume imaginary part
> = 0). Antonio's functions handle complex & non-complex AFAICS
>
>> , and will also expose the derived band somehow (can I
>> benefit from the CreateDerivedBand API,
> If you write it, yes... (it doesn't exist yet if I wasn't clear)
>
> If you follow the VRT pixel function road, that's a matter of building a on-
> the-fly VRT dataset.
>
> http://www.gdal.org/gdal_vrttut.html#gdal_vrttut_creation could possibly be
> used. But AFAICS setting the pixel function name isn't supported. Could be
> done by implementing CPLErr VRTSourcedRasterBand::SetMetadataItem( const char
> *pszName, const char *pszValue, const char *pszDomain ), similarly to what
> VRTSourcedRasterBand::SetMetadataItem(...) handles for vrt sources.
>
>> and the pixel functions from Antonio ?).
> Would make sense to internalize this code (uses X/MIT license)
>
>> Also, "original_datasetname" can be a subdataset itself.
> Yes, that shouldn't matter to the driver. It is just a string that it provides
> to GDALOpen() (some care must be taken when parsing
> DERIVED_SUBDATASET:algname:original_datasetname, since original_datasetname
> can have columns too)
>
>> There should also be a mechanism that will report the existing
>> DERIVED_SUBDATASETS upon query (this I did not get yet).
> Implement in GDALDataset, the following method (from GDALMajorObject)
>
> virtual char **GetMetadata( const char * pszDomain = "" );
>
> char **GDALDataset::GetMetadata( const char * pszDomain = "" )
> {
> if( pszDomain != NULL && EQUAL(pszDomain, "DERIVED_SUBDATASET") )
> {
> // some stuff. Be careful that ownership of returned char** belongs
> // to the dataset object
> }
> else
> return GDALMajorObject::GetMetadata(pszDomain);
> }
>
> As well as :
>
> virtual char **GetMetadataDomainList(); to report DERIVED_SUBDATASET when it
> makes sense (contrary to GetMetadata(), ownership of the returned char**
> belongs to the caller)
>
> Even
>


--
Julien MICHEL
CNES - DCT/SI/AP - BPI 1219
18, avenue Edouard Belin
31401 Toulouse Cedex 09 - France
Tel: +33 561 282 894 - Fax: +33 561 283 109

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

Re: Amplitude virtual bands for complex datasets ?

Julien Michel-2
PS: any comments welcome, I am not really familiar with gdal coding
style & arch.

Julien

Le 04/07/2016 à 14:30, Julien Michel a écrit :

> Hi,
>
> Started the work here :
>
> https://github.com/OSGeo/gdal/compare/trunk...jmichel-otb:enhance-complex-datasets 
>
>
> For now, recognizes the following syntax :
>
> gdal_translate -srcwin 1000 1000 1000 1000
> DERIVED_SUBDATASET:COMPLEX_AMPLITUDE:s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff
> ~/tmp/test.tif
>
> I need to see about exposing those derived datasets properly.
>
> Antonio, I copied one of your pixel function for fast prototyping, I
> hope you do not mind. Of course upon completion of this work we should
> merge properly all your pixel functions with proper
> credits/copyright/licence. Thing is I did not know where to put them
> within gdal.
>
> Regards,
>
> Julien
>
> Le 04/07/2016 à 11:06, Even Rouault a écrit :
>> Hi,
>>> To check if I understood well :
>>>
>>> I will create a driver that will recognize the
>>> "DERIVED_SUBDATASET:Amplitude:original_datasetname" syntax. This driver
>>> needs to now if "original_datasetname" is of complex type (to report it
>>> can open it)
>> That or we could also possibly accept non-complex types (assume
>> imaginary part
>> = 0). Antonio's functions handle complex & non-complex AFAICS
>>
>>> , and will also expose the derived band somehow (can I
>>> benefit from the CreateDerivedBand API,
>> If you write it, yes... (it doesn't exist yet if I wasn't clear)
>>
>> If you follow the VRT pixel function road, that's a matter of
>> building a on-
>> the-fly VRT dataset.
>>
>> http://www.gdal.org/gdal_vrttut.html#gdal_vrttut_creation could
>> possibly be
>> used. But AFAICS setting the pixel function name isn't supported.
>> Could be
>> done by implementing CPLErr VRTSourcedRasterBand::SetMetadataItem(
>> const char
>> *pszName, const char *pszValue, const char *pszDomain ), similarly to
>> what
>> VRTSourcedRasterBand::SetMetadataItem(...) handles for vrt sources.
>>
>>> and the pixel functions from Antonio ?).
>> Would make sense to internalize this code (uses X/MIT license)
>>
>>> Also, "original_datasetname" can be a subdataset itself.
>> Yes, that shouldn't matter to the driver. It is just a string that it
>> provides
>> to GDALOpen() (some care must be taken when parsing
>> DERIVED_SUBDATASET:algname:original_datasetname, since
>> original_datasetname
>> can have columns too)
>>
>>> There should also be a mechanism that will report the existing
>>> DERIVED_SUBDATASETS upon query (this I did not get yet).
>> Implement in GDALDataset, the following method (from GDALMajorObject)
>>
>> virtual char **GetMetadata( const char * pszDomain = "" );
>>
>> char **GDALDataset::GetMetadata( const char * pszDomain = "" )
>> {
>>     if( pszDomain != NULL && EQUAL(pszDomain, "DERIVED_SUBDATASET") )
>>     {
>>         // some stuff. Be careful that ownership of returned char**
>> belongs
>>         // to the dataset object
>>     }
>>     else
>>         return GDALMajorObject::GetMetadata(pszDomain);
>> }
>>
>> As well as :
>>
>> virtual char **GetMetadataDomainList(); to report DERIVED_SUBDATASET
>> when it
>> makes sense (contrary to GetMetadata(), ownership of the returned char**
>> belongs to the caller)
>>
>> Even
>>
>
>


--
Julien MICHEL
CNES - DCT/SI/AP - BPI 1219
18, avenue Edouard Belin
31401 Toulouse Cedex 09 - France
Tel: +33 561 282 894 - Fax: +33 561 283 109

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

Re: Amplitude virtual bands for complex datasets ?

Even Rouault-2
In reply to this post by Julien Michel-2
Le lundi 04 juillet 2016 14:30:26, Julien Michel a écrit :

> Hi,
>
> Started the work here :
>
> https://github.com/OSGeo/gdal/compare/trunk...jmichel-otb:enhance-complex-d
> atasets
>
> For now, recognizes the following syntax :
>
> gdal_translate -srcwin 1000 1000 1000 1000
> DERIVED_SUBDATASET:COMPLEX_AMPLITUDE:s1a-s6-slc-vv-20150619t195043-20150619
> t195101-006447-00887d-001.tiff ~/tmp/test.tif
>
> I need to see about exposing those derived datasets properly.
>
> Antonio, I copied one of your pixel function for fast prototyping, I
> hope you do not mind. Of course upon completion of this work we should
> merge properly all your pixel functions with proper
> credits/copyright/licence. Thing is I did not know where to put them
> within gdal.

frmts/vrt seem to be the appropriate place, and make the VRTDriver register
them.

>
> Regards,
>
> Julien
>
> Le 04/07/2016 à 11:06, Even Rouault a écrit :
> > Hi,
> >
> >> To check if I understood well :
> >>
> >> I will create a driver that will recognize the
> >> "DERIVED_SUBDATASET:Amplitude:original_datasetname" syntax. This driver
> >> needs to now if "original_datasetname" is of complex type (to report it
> >> can open it)
> >
> > That or we could also possibly accept non-complex types (assume imaginary
> > part = 0). Antonio's functions handle complex & non-complex AFAICS
> >
> >> , and will also expose the derived band somehow (can I
> >> benefit from the CreateDerivedBand API,
> >
> > If you write it, yes... (it doesn't exist yet if I wasn't clear)
> >
> > If you follow the VRT pixel function road, that's a matter of building a
> > on- the-fly VRT dataset.
> >
> > http://www.gdal.org/gdal_vrttut.html#gdal_vrttut_creation could possibly
> > be used. But AFAICS setting the pixel function name isn't supported.
> > Could be done by implementing CPLErr
> > VRTSourcedRasterBand::SetMetadataItem( const char *pszName, const char
> > *pszValue, const char *pszDomain ), similarly to what
> > VRTSourcedRasterBand::SetMetadataItem(...) handles for vrt sources.
> >
> >> and the pixel functions from Antonio ?).
> >
> > Would make sense to internalize this code (uses X/MIT license)
> >
> >> Also, "original_datasetname" can be a subdataset itself.
> >
> > Yes, that shouldn't matter to the driver. It is just a string that it
> > provides to GDALOpen() (some care must be taken when parsing
> > DERIVED_SUBDATASET:algname:original_datasetname, since
> > original_datasetname can have columns too)
> >
> >> There should also be a mechanism that will report the existing
> >> DERIVED_SUBDATASETS upon query (this I did not get yet).
> >
> > Implement in GDALDataset, the following method (from GDALMajorObject)
> >
> > virtual char **GetMetadata( const char * pszDomain = "" );
> >
> > char **GDALDataset::GetMetadata( const char * pszDomain = "" )
> > {
> >
> > if( pszDomain != NULL && EQUAL(pszDomain, "DERIVED_SUBDATASET") )
> > {
> >
> > // some stuff. Be careful that ownership of returned char** belongs
> > // to the dataset object
> >
> > }
> > else
> >
> > return GDALMajorObject::GetMetadata(pszDomain);
> >
> > }
> >
> > As well as :
> >
> > virtual char **GetMetadataDomainList(); to report DERIVED_SUBDATASET when
> > it makes sense (contrary to GetMetadata(), ownership of the returned
> > char** belongs to the caller)
> >
> > Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Amplitude virtual bands for complex datasets ?

Antonio Valentino
In reply to this post by Julien Michel-2
Hi Julien,

Il 04/07/2016 14:30, Julien Michel ha scritto:

> Hi,
>
> Started the work here :
>
> https://github.com/OSGeo/gdal/compare/trunk...jmichel-otb:enhance-complex-datasets
>
>
> For now, recognizes the following syntax :
>
> gdal_translate -srcwin 1000 1000 1000 1000
> DERIVED_SUBDATASET:COMPLEX_AMPLITUDE:s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff
> ~/tmp/test.tif
>
> I need to see about exposing those derived datasets properly.
>
> Antonio, I copied one of your pixel function for fast prototyping, I
> hope you do not mind. Of course upon completion of this work we should
absolutely, I'm happy if it helps

> merge properly all your pixel functions with proper
> credits/copyright/licence. Thing is I did not know where to put them
> within gdal.


cheers

--
Antonio Valentino


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

signature.asc (188 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Amplitude virtual bands for complex datasets ?

Bugbuster
In reply to this post by Even Rouault-2
Dear Even,

Sorry for my late answer, I was onto something different yesterday.

I agree with you, and I did change all function declarations in order to "export" them in the DLL.
Thus, all functions in both pixfunplugin.c and pixelfunctions.c are now declared as "C" functions with the CPL_DLL export tag.
I checked with Dependency Walker that all functions are indeed exported correctly.
What is really weird is that all I/F variables in RealPixFunc() function (for example) seem to be correct but the 1sr one papoSources.

Philippe
Reply | Threaded
Open this post in threaded view
|

Re: Amplitude virtual bands for complex datasets ?

Julien Michel-2
In reply to this post by Even Rouault-2
Hi all,

I hacked the SetBand() method so as to report the new medata in
DERIVED_SUBDATASETS mdd. Now if I gdal info a S1 dataset I get :

$ gdalinfo -nogcp -mdd DERIVED_SUBDATASETS
s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff
Driver: GTiff/GeoTIFF
Files: s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff
s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff.aux.xml
Size is 17663, 31106
[...]
Metadata (DERIVED_SUBDATASETS):
   DERIVED_SUBDATASET_1_DESC=Complex amplitude of bands from
/media/michelj/HUGE/LargeInput/SENTINEL1/S1A_S6_SLC__1SSV_20150619T195043/measurement/s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff
DERIVED_SUBDATASET_1_NAME=DERIVED_SUBDATASET:COMPLEX_AMPLITUDE:/media/michelj/HUGE/LargeInput/SENTINEL1/S1A_S6_SLC__1SSV_20150619T195043/measurement/s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff
[...]
Band 1 Block=17663x1 Type=CInt16, ColorInterp=Gray

Even, I think the AddBand() method is obviously not the right place to
do that, but I do not get how tweak the GetMetadata() /
GetMetadataDomainList(). As far as I can tell from the code, datasets
never directly declare new domains by hacking GetMetadataDomainList().
Calling SetMetadataItem() with a new domain is enough to create it (that
is what I did). But I surely miss something here.

Also, can we discuss syntax (COMPLEX_AMPLITUDE is not very meaningful)
and behaviour (for now all bands are considered even if they are not
complex).

Last, it does not work for real subdatasets. For instance, for S2 data :

$ gdalinfo -mdd DERIVED_SUBDATASETS
SENTINEL2_L1C:S2A_OPER_MTD_SAFL1C_PDMC_20150519T113602_R070_V20130707T172156_20130707T172156.xml:20m:EPSG_32615

Driver: SENTINEL2/Sentinel 2
[...]
Metadata (DERIVED_SUBDATASETS):
   DERIVED_SUBDATASET_1_DESC=Complex amplitude of bands from
   DERIVED_SUBDATASET_1_NAME=DERIVED_SUBDATASET:COMPLEX_AMPLITUDE:
[...]

But calling
$ gdalinfo
DERIVED_SUBDATASET:COMPLEX_AMPLITUDE:SENTINEL2_L1C:S2A_OPER_MTD_SAFL1C_PDMC_20150519T113602_R070_V20130707T172156_20130707T172156.xml:20m:EPSG_32615

seems to work.

I will try to do the pixel function proper integration on a separate
branch, so as to avoid mixing different purposes.

Regards,

Julien

Le 04/07/2016 à 15:04, Even Rouault a écrit :

> Le lundi 04 juillet 2016 14:30:26, Julien Michel a écrit :
>> Hi,
>>
>> Started the work here :
>>
>> https://github.com/OSGeo/gdal/compare/trunk...jmichel-otb:enhance-complex-d
>> atasets
>>
>> For now, recognizes the following syntax :
>>
>> gdal_translate -srcwin 1000 1000 1000 1000
>> DERIVED_SUBDATASET:COMPLEX_AMPLITUDE:s1a-s6-slc-vv-20150619t195043-20150619
>> t195101-006447-00887d-001.tiff ~/tmp/test.tif
>>
>> I need to see about exposing those derived datasets properly.
>>
>> Antonio, I copied one of your pixel function for fast prototyping, I
>> hope you do not mind. Of course upon completion of this work we should
>> merge properly all your pixel functions with proper
>> credits/copyright/licence. Thing is I did not know where to put them
>> within gdal.
> frmts/vrt seem to be the appropriate place, and make the VRTDriver register
> them.
>
>> Regards,
>>
>> Julien
>>
>> Le 04/07/2016 à 11:06, Even Rouault a écrit :
>>> Hi,
>>>
>>>> To check if I understood well :
>>>>
>>>> I will create a driver that will recognize the
>>>> "DERIVED_SUBDATASET:Amplitude:original_datasetname" syntax. This driver
>>>> needs to now if "original_datasetname" is of complex type (to report it
>>>> can open it)
>>> That or we could also possibly accept non-complex types (assume imaginary
>>> part = 0). Antonio's functions handle complex & non-complex AFAICS
>>>
>>>> , and will also expose the derived band somehow (can I
>>>> benefit from the CreateDerivedBand API,
>>> If you write it, yes... (it doesn't exist yet if I wasn't clear)
>>>
>>> If you follow the VRT pixel function road, that's a matter of building a
>>> on- the-fly VRT dataset.
>>>
>>> http://www.gdal.org/gdal_vrttut.html#gdal_vrttut_creation could possibly
>>> be used. But AFAICS setting the pixel function name isn't supported.
>>> Could be done by implementing CPLErr
>>> VRTSourcedRasterBand::SetMetadataItem( const char *pszName, const char
>>> *pszValue, const char *pszDomain ), similarly to what
>>> VRTSourcedRasterBand::SetMetadataItem(...) handles for vrt sources.
>>>
>>>> and the pixel functions from Antonio ?).
>>> Would make sense to internalize this code (uses X/MIT license)
>>>
>>>> Also, "original_datasetname" can be a subdataset itself.
>>> Yes, that shouldn't matter to the driver. It is just a string that it
>>> provides to GDALOpen() (some care must be taken when parsing
>>> DERIVED_SUBDATASET:algname:original_datasetname, since
>>> original_datasetname can have columns too)
>>>
>>>> There should also be a mechanism that will report the existing
>>>> DERIVED_SUBDATASETS upon query (this I did not get yet).
>>> Implement in GDALDataset, the following method (from GDALMajorObject)
>>>
>>> virtual char **GetMetadata( const char * pszDomain = "" );
>>>
>>> char **GDALDataset::GetMetadata( const char * pszDomain = "" )
>>> {
>>>
>>> if( pszDomain != NULL && EQUAL(pszDomain, "DERIVED_SUBDATASET") )
>>> {
>>>
>>> // some stuff. Be careful that ownership of returned char** belongs
>>> // to the dataset object
>>>
>>> }
>>> else
>>>
>>> return GDALMajorObject::GetMetadata(pszDomain);
>>>
>>> }
>>>
>>> As well as :
>>>
>>> virtual char **GetMetadataDomainList(); to report DERIVED_SUBDATASET when
>>> it makes sense (contrary to GetMetadata(), ownership of the returned
>>> char** belongs to the caller)
>>>
>>> Even


--
Julien MICHEL
CNES - DCT/SI/AP - BPI 1219
18, avenue Edouard Belin
31401 Toulouse Cedex 09 - France
Tel: +33 561 282 894 - Fax: +33 561 283 109

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

Re: Amplitude virtual bands for complex datasets ?

Even Rouault-2
>
> Even, I think the AddBand() method is obviously not the right place to
> do that, but I do not get how tweak the GetMetadata() /
> GetMetadataDomainList().

Hum I'm not sure how to better explain. I almost coded it ;-)

> As far as I can tell from the code, datasets
> never directly declare new domains by hacking GetMetadataDomainList().
> Calling SetMetadataItem() with a new domain is enough to create it (that
> is what I did). But I surely miss something here.

We definitely don't want to call SetMetadata() / SetMetadataItem() in generic
code, as it could have a lot of undesired side effects, like pushing metadata
to PAM (.aux.xml files) or to the own metadata storage of the driver. IMHO, the
derived subdataset reporting must rather be done in a "passive" way with
GetMetadata().

Hum, after all, maybe SetBand() could still do it, but we would need to make
sure to explicitly call GDALDataset::SetMetadataItem(), and not the potential
overridden SetMetadataItem() method of a driver.

>
> Also, can we discuss syntax (COMPLEX_AMPLITUDE is not very meaningful)
> and behaviour (for now all bands are considered even if they are not
> complex).
>
> Last, it does not work for real subdatasets. For instance, for S2 data :
>
> $ gdalinfo -mdd DERIVED_SUBDATASETS
> SENTINEL2_L1C:S2A_OPER_MTD_SAFL1C_PDMC_20150519T113602_R070_V20130707T17215
> 6_20130707T172156.xml:20m:EPSG_32615
>
> Driver: SENTINEL2/Sentinel 2
> [...]
> Metadata (DERIVED_SUBDATASETS):
>    DERIVED_SUBDATASET_1_DESC=Complex amplitude of bands from
>    DERIVED_SUBDATASET_1_NAME=DERIVED_SUBDATASET:COMPLEX_AMPLITUDE:
> [...]

Probably because at the time of SetBand(), the dataset name hasn't been set
yet on the dataset object. It is generally done at the very end of the Open()
method, and that's the case of the Sentinel2 driver, and most all other
drivers that are generally copied&pasted from the same skeleton. So the
SetBand() way seems to be a dead end, unless we edit all drivers (probably not
all drivers, just in drivers that can have complex data types) to initialize
their description earlier.


Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Amplitude virtual bands for complex datasets ?

Julien Michel-2
Le 06/07/2016 à 10:26, Even Rouault a écrit :
>> Even, I think the AddBand() method is obviously not the right place to
>> do that, but I do not get how tweak the GetMetadata() /
>> GetMetadataDomainList().
> Hum I'm not sure how to better explain. I almost coded it ;-)
Sorry, it is my first time coding into gdal, I need to wrap my head
around it. I will try again ;)

>
>> As far as I can tell from the code, datasets
>> never directly declare new domains by hacking GetMetadataDomainList().
>> Calling SetMetadataItem() with a new domain is enough to create it (that
>> is what I did). But I surely miss something here.
> We definitely don't want to call SetMetadata() / SetMetadataItem() in generic
> code, as it could have a lot of undesired side effects, like pushing metadata
> to PAM (.aux.xml files) or to the own metadata storage of the driver. IMHO, the
> derived subdataset reporting must rather be done in a "passive" way with
> GetMetadata().
Pushing metadata to PAM actually happens with the current code. I had a
hard time figuring out why my changes to the derived datasets name did
not reflect in gdalinfo, before I remember this little xml file ;) I
understand this is not a desired behaviour.
>
> Probably because at the time of SetBand(), the dataset name hasn't been set
> yet on the dataset object. It is generally done at the very end of the Open()
> method, and that's the case of the Sentinel2 driver, and most all other
> drivers that are generally copied&pasted from the same skeleton. So the
> SetBand() way seems to be a dead end, unless we edit all drivers (probably not
> all drivers, just in drivers that can have complex data types) to initialize
> their description earlier.
>
Ok I will try the other way, and keep you posted.

Regards,

Julien

--
Julien MICHEL
CNES - DCT/SI/AP

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