i.atcorr: Sentinel-2 support broken on some systems

[sbl]

It seems that Sentinel-2 functions in i.atcorr are effectively broken in some release packages and on some systems.

The following installations have been reported to only return NULL values in output:

• Ubuntu 16.04 with GRASS 7.4.0RC1 from UbuntuGIS-experimental
• UBUNTU 14.04 with GRASS 7.4.0RC1 self compiled with gcc 5.4.1 (and older)
• UBUNTU 14.04 with GRASS 7.5 self compiled with gcc 4.8
• Windows 8.1 with GRASS 7.4.0RC1 and 7.5 daily from OSGeo4W

However, the following systems produce reasonable output:

• UBUNTU 16.04.3 LTS with GRASS 7.4 and 7.5 self compiled with gcc 5.4.0 20160609 Ubuntu 5.4.0-6ubuntu1~16.04.5
• Fedora 27 with GRASS 7.4.0RC1 and 7.5 self compiled with gcc 7.2.1 20170915

Test case below, test data attached to the ticket.

r.in.gdal input=dem.tif output=dem –o -–v --o
r.in.gdal input=S2A_OPER_PRD_MSIL1C_PDMC_20160907T044118_R008_V20160905T104022_20160905T104245.B08.tif output=B08 -o -–v --o
g.region -p raster=B08 align=B08
i.atcorr -r input=B08 elevation=dem range=0,10000 output=test_atcorr rescale=0,10000 parameters=p6s.txt -–v --o

relevant compiler flags are:

-g -Wall -Wextra -Wpedantic -Wshadow -Wno-sign-compare -fno-common -fexceptions -Werror=implicit-function-declaration -Wp,-D_FORTIFY_SOURCE=2 -O3 -fno-fast-math

-fno-fast-math might be important

Finally, running i.atcorr through valgrind (on Fedora) gives lots of the following warnings, that according to MarkusM should be investigated:

==14080== Conditional jump or move depends on uninitialised value(s)
==14080==    at 0x42B18D: os(float, float, float, float, float, float (&) [51][49], Gauss&, Altitude const&, GeomCond const&) (computations.cpp:831)
==14080==    by 0x42E9E3: atmref(float, float, float, float, float, OpticalAtmosProperties&, Gauss&, GeomCond const&, AerosolModel const&, Altitude const&) (computations.cpp:1408)
==14080==    by 0x42FCBC: discom(GeomCond const&, AtmosModel const&, AerosolModel const&, AerosolConcentration const&, Altitude const&, IWave const&) (computations.cpp:1654)
==14080==    by 0x40C032: init_6S(char*) (6s.cpp:100)
==14080==    by 0x405A3C: main (main.cpp:618)
...
later on

==14080== Conditional jump or move depends on uninitialised value(s)
==14080==    at 0x42DED8: iso(float, float, float, float, float, float (&) [3], Gauss&, Altitude const&) (computations.cpp:1262)
==14080==    by 0x42F011: scatra(float, float, float, float, float, OpticalAtmosProperties&, Gauss&, GeomCond const&, Altitude const&) (computations.cpp:1578)
==14080==    by 0x42FD08: discom(GeomCond const&, AtmosModel const&, AerosolModel const&, AerosolConcentration const&, Altitude const&, IWave const&) (computations.cpp:1659)
==14080==    by 0x406882: pre_compute_h(float) (6s.cpp:148)
==14080==    by 0x406372: process_raster (main.cpp:369)
==14080==    by 0x406372: main (main.cpp:632)
...

For ML discussion see: ​https://lists.osgeo.org/pipermail/grass-user/2017-December/077545.html

[sbl]

6s parameters (UNIX line endings)

[sbl]

DEM

[sbl]

Sentinel-2 Band 8

[mlennert]

I'm on Debian testing with gcc (Debian 7.2.0-17) 7.2.1 20171205.

Using fresh trunk with the following compiler settings:

CFLAGS="-g -Wall -fopenmp -lgomp -Ofast -fno-fast-math
-march=core-avx-i -fno-common -fexceptions $PENTIUM64"
LDFLAGS="-Wl,--no-undefined -fopenmp -lgomp"
CXXFLAGS="-g -Ofast"

and applying the exact commands given in this ticket, I get as result:

r.info -r test_atcorr
min=NULL
max=NULL

Several questions:

• In the discussion in the mailing list, some examples were with i.atcorr -r, others without. Do we agree that in this case it should be with -r as the data contains reflectance information ?
• What is the exact meaning and impact of the range setting ?
• Stefan, any special reason why you import your data with '-o' ? Why not create an EPSG 32633 location and import the data there ?

[mmetz]

Replying to mlennert:

I'm on Debian testing with gcc (Debian 7.2.0-17) 7.2.1 20171205.

Using fresh trunk with the following compiler settings:

CFLAGS="-g -Wall -fopenmp -lgomp -Ofast -fno-fast-math
-march=core-avx-i -fno-common -fexceptions $PENTIUM64"
LDFLAGS="-Wl,--no-undefined -fopenmp -lgomp"
CXXFLAGS="-g -Ofast"

-Ofast is dangerous. I would regard any optimization larger than -O3 as experimental. Some code only runs properly with max -O1.

[sbl]

Replying to mlennert:

Several questions:

• In the discussion in the mailing list, some examples were with i.atcorr -r, others without. Do we agree that in this case it should be with -r as the data contains reflectance information ?

Yes, I would say so.

• Stefan, any special reason why you import your data with '-o' ? Why not create an EPSG 32633 location and import the data there ?

Laziness ;-): We usually treat 32633 and 25833 as identical (max difference in Norway is reported to be ~40cm). So, I just did not want to create a new location for testing with Sajids data...

[mlennert]

Replying to mmetz:

Replying to mlennert:

I'm on Debian testing with gcc (Debian 7.2.0-17) 7.2.1 20171205.

Using fresh trunk with the following compiler settings:

CFLAGS="-g -Wall -fopenmp -lgomp -Ofast -fno-fast-math
-march=core-avx-i -fno-common -fexceptions $PENTIUM64"
LDFLAGS="-Wl,--no-undefined -fopenmp -lgomp"
CXXFLAGS="-g -Ofast"

-Ofast is dangerous. I would regard any optimization larger than -O3 as experimental. Some code only runs properly with max -O1.

I just recompiled with all compiler flags commented out. Same result, i.e. all NULLs. Now I have to go into a meeting and cannot test further.

[mlennert]

Using

CFLAGS=-g -Wall -Wextra -Wpedantic -Wshadow -Wno-sign-compare -fno-common -fexceptions -Werror=implicit-function-declaration -Wp,-D_FORTIFY_SOURCE=2 -O3 -fno-fast-math

I again get

r.info -r test_atcorr
min=NULL
max=NULL

So not sure it's the compiler flags...

[mmetz]

Replying to mlennert:

Using

CFLAGS=-g -Wall -Wextra -Wpedantic -Wshadow -Wno-sign-compare -fno-common -fexceptions -Werror=implicit-function-declaration -Wp,-D_FORTIFY_SOURCE=2 -O3 -fno-fast-math

Are you still using

CXXFLAGS="-g -Ofast"

?

i.atcorr is a C++ module

I again get

r.info -r test_atcorr
min=NULL
max=NULL

So not sure it's the compiler flags...

I can confirm on Debian testing, all NULL...

[mlennert]

Replying to mmetz:

Replying to mlennert:

Using

CFLAGS=-g -Wall -Wextra -Wpedantic -Wshadow -Wno-sign-compare -fno-common -fexceptions -Werror=implicit-function-declaration -Wp,-D_FORTIFY_SOURCE=2 -O3 -fno-fast-math

Are you still using

CXXFLAGS="-g -Ofast"

?

No. I left all ldflags and cxxflags commented, and I don't think -Ofast is default.

i.atcorr is a C++ module

Good to know. So do CFLAGS have any impact ?

I again get

r.info -r test_atcorr
min=NULL
max=NULL

So not sure it's the compiler flags...

I can confirm on Debian testing, all NULL...

Good. This means that it is something systematic which we should be able to find. Stefan, you tried different versions. Are there any older versions which work for you ?

[mmetz]

Replying to mlennert:

[...] This means that it is something systematic which we should be able to find. Stefan, you tried different versions. Are there any older versions which work for you ?

It seems that the 6s code in i.atcorr is numerically unstable: some internal variables can get very close to zero. If these close-to-zero values are snapped to zero, nan and inf results can appear. This is why the result can be all NULL, or not.

The symptom can be cured by using double precision floating point variables throughout. Currently, i.atcorr uses single precision floating point variables. When using double precision floating point variables, Debian testing produces reasonable results. However, these results differ from Fedora 27, also using double precision floating point variables throughout. I would say that both results (Debian testing and Fedora 27) are slightly wrong because of the numerical instability of the 6s code in combination with the sentinel-2 parameters.

I suggest to review the parameters for sentinel-2 as used by i.atcorr.

I also suggest to review i.atcorr. The code has been translated from Fortran to C++ and might need some manual adjustments.

[sbl]

Replying to mmetz:

The symptom can be cured by using double precision floating point variables throughout. Currently, i.atcorr uses single precision floating point variables. When using double precision floating point variables, Debian testing produces reasonable results. However, these results differ from Fedora 27, also using double precision floating point variables throughout. I would say that both results (Debian testing and Fedora 27) are slightly wrong because of the numerical instability of the 6s code in combination with the sentinel-2 parameters.

Would you be able to provide a patch with those changes? I would not mind slightly wrong results, just I am able to proceed. I can update the computation once this is properly fixed...

I suggest to review the parameters for sentinel-2 as used by i.atcorr. I also suggest to review i.atcorr. The code has been translated from Fortran to C++ and might need some manual adjustments.

I will happily test any fixes on various systems!

[mmetz]

patch to use double precision fp

[mmetz]

Replying to sbl:

Replying to mmetz:

The symptom can be cured by using double precision floating point variables throughout. Currently, i.atcorr uses single precision floating point variables. When using double precision floating point variables, Debian testing produces reasonable results. However, these results differ from Fedora 27, also using double precision floating point variables throughout. I would say that both results (Debian testing and Fedora 27) are slightly wrong because of the numerical instability of the 6s code in combination with the sentinel-2 parameters.

Would you be able to provide a patch with those changes? I would not mind slightly wrong results, just I am able to proceed. I can update the computation once this is properly fixed...

Patch i_atcorr_double.patch attached. i.atcorr is identical in 7.4 and 7.5, therefore the patch works with both versions.

[mlennert]

Replying to mmetz:

Replying to sbl:

Replying to mmetz:

The symptom can be cured by using double precision floating point variables throughout. Currently, i.atcorr uses single precision floating point variables. When using double precision floating point variables, Debian testing produces reasonable results. However, these results differ from Fedora 27, also using double precision floating point variables throughout. I would say that both results (Debian testing and Fedora 27) are slightly wrong because of the numerical instability of the 6s code in combination with the sentinel-2 parameters.

Would you be able to provide a patch with those changes? I would not mind slightly wrong results, just I am able to proceed. I can update the computation once this is properly fixed...

Patch i_atcorr_double.patch attached. i.atcorr is identical in 7.4 and 7.5, therefore the patch works with both versions.

Thanks for investigating this !

I now get what looks like more reasonable results:

r.info -r test_atcorr
min=0.3878489
max=8513.962

compared to the original before correction:

r.info -r B08
min=83
max=7110

Interestingly, i.atcorr also runs much faster with the patch.

[sbl]

Yes, thanks so much Markus! With the patch i.atcorr gives now also numbers on Ubuntu 14.04:

min = 0  max = 8514

So results are identical to Moritz`s.

I noticed that although "rescale" option was set to be 0-10000, maximum output was 8514...

With e.g. Landsat paramters, i.atcorr scales to the full range provided in rescale...

And yes, it is much faster for Sentinel-2 now...

[mmetz]

Replying to sbl:

Yes, thanks so much Markus! With the patch i.atcorr gives now also numbers on Ubuntu 14.04:

min = 0  max = 8514

So results are identical to Moritz`s.

I noticed that although "rescale" option was set to be 0-10000, maximum output was 8514...

With

i.atcorr -r input=B08 elevation=dem range=0,10000 output=test_atcorr rescale=0,10000 parameters=p6s.txt

I get

min=8.736016e-10
max=10000

on Fedora 27

and

min=8.708314e-10
max=10000

on Debian testing (Buster)

So there is some non-NULL output, but it differs between systems. Not so nice.

[sbl]

Replying to mmetz:

So there is some non-NULL output, but it differs between systems. Not so nice.

Sorry, we should have reported the command we used. If I run exactly your case:

i.atcorr -r input=B08 elevation=dem range=0,10000 output=test_atcorr rescale=0,10000 parameters=p6s.txt

I get the same output as you:

r.info -r test_atcorr
min=8.708314e-10
max=10000

So, it seems to be essential to use both range and rescale option if one wants the data to be exactly in a given range of values... Maybe something to clarify in the manual or set as a parser rule / warning in the module?

[mlennert]

Replying to sbl:

Yes, thanks so much Markus! With the patch i.atcorr gives now also numbers on Ubuntu 14.04:

min = 0  max = 8514

So results are identical to Moritz`s.

I noticed that although "rescale" option was set to be 0-10000, maximum output was 8514...

This depends on your input data and your range parameter:

r.info -r B08
min=83
max=7110

i.atcorr -r input=B08 elevation=dem range=83,7110 output=test_atcorr rescale=0,10000 parameters=BROL/STEFAN/p6s.txt --v --o

r.info -r test_atcorr
min=8.732422e-10
max=10000

i.e. if you tell it to scale the input data range to 0,10000 it will do so, but if you tell it that the input range is 1,10000 then it will scale this range to 0,10000.

Markus, could it be that you are using the original larger input image from the ML and not the smaller clip attached to this ticket ?

[sbl]

Replying to mlennert:

Markus, could it be that you are using the original larger input image from the ML and not the smaller clip attached to this ticket ?

Yes, thanks Moritz! That was it in my case! Now:

r.info -r B08
min=83
max=7110

i.atcorr -r input=B08 elevation=dem range=83,7110 output=test_atcorr rescale=0,10000 parameters=p6s.txt --o
Atmospheric correction...
 100%
Atmospheric correction complete.
r.info -r test_atcorr
min=8.732422e-10
max=10000

and

i.atcorr -r input=B08 elevation=dem output=test_atcorr rescale=0,10000 parameters=p6s.txt --o
Atmospheric correction...
 100%
Atmospheric correction complete.
r.info -r test_atcorr
min=3875.138
max=10000

[mmetz]

Replying to mlennert:

Replying to sbl:

Yes, thanks so much Markus! With the patch i.atcorr gives now also numbers on Ubuntu 14.04:

min = 0  max = 8514

So results are identical to Moritz`s.

I noticed that although "rescale" option was set to be 0-10000, maximum output was 8514...

This depends on your input data and your range parameter:

r.info -r B08
min=83
max=7110

i.atcorr -r input=B08 elevation=dem range=83,7110 output=test_atcorr rescale=0,10000 parameters=BROL/STEFAN/p6s.txt --v --o

r.info -r test_atcorr
min=8.732422e-10
max=10000

i.e. if you tell it to scale the input data range to 0,10000 it will do so, but if you tell it that the input range is 1,10000 then it will scale this range to 0,10000.

Markus, could it be that you are using the original larger input image from the ML and not the smaller clip attached to this ticket ?

Yes, right. Now with the smaller clip using the command in the description of this ticket

i.atcorr -r input=B08 elevation=dem range=0,10000 output=test_atcorr rescale=0,10000 parameters=p6s.txt --v --o

I get

min=0.04157807
max=8513.816

you get

min=0.3878489
max=8513.962

i.atcorr rescales the input to 0,1 using the range option, then does the transformation and rescales the result from 0,1 to the rescale option.

Be aware that if you don't specify range and rescale, default 0,255 will be used which will produce strange results if the input range is much larger than 0,255. Conversely, if the actual input range is smaller than provided with the range option, the output will also be only a subset of the rescale range.

[sbl]

Would it be an option to apply the patch to trunk / 7.4.0RC2 ?

Or do you prefer "no output" (on Debian based systems) over "silent and possibly slightly wrong" output? I can see that the patch might affect results for other sensors, but it is mainly a precision issue, isn`t it?

What tests would be needed in order to get a production ready version of Sentinel-2 in i.atcorr?

Do we need to test if other sensors are affected significantly by this patch?

Any idea how to find out what is causing the difference between systems? Does not seem to be compiler flags...

[mmetz]

patch to fix reading spectral responses

[mmetz]

Replying to sbl:

Would it be an option to apply the patch to trunk / 7.4.0RC2 ?

No, the patch was more meant as a diagnostic tool.

The Sentinel-2 parameters are somewhat unusual with many zeroes in the spectral responses. For othjer sensors, the leading and trailing zeroes have been clipped off. The new patch i_atcoor_iwave.patch fixes these leading and trailing zeroes on initialization. This patch affects all sensors, not only Sentinel-2. There have been regular reports that the output of i.atcorr is all NULL, also for other sensors. This patch could fix this, and the previous patch converting float to double is no longer needed.

[sbl]

Thanks Markus for the swift reply! The patch works fine:

g.version -gr
version=7.4.0svn
date=2018
revision=r72045M
build_date=2018-01-08
build_platform=x86_64-pc-linux-gnu
build_off_t_size=8
libgis_revision=70829 
libgis_date="2017-04-04 09:43:02 +0200 (Tue, 04 Apr 2017) "

r.info -r B08
min=83
max=7110

g.region raster=B08 align=B08 -p
projection: 1 (UTM)
zone:       33
datum:      etrs89
ellipsoid:  grs80
north:      6649695
south:      6647685
west:       260495
east:       262505
nsres:      10
ewres:      10
rows:       201
cols:       201
cells:      40401

i.atcorr -r input=B08 elevation=dem range=0,10000 output=test_atcorr rescale=0,10000 parameters=~/p6s.txt --v --o
(...)
Atmospheric correction...
 100%
Atmospheric correction complete.

r.info -r test_atcorr
min=0.5415085
max=8514.146

[mmetz]

Replying to sbl:

With the patch applied to clean and completely fresh svn checkout I get the following results (on UBUNTU 14.04):

[...]

i.atcorr -r input=B08 elevation=dem range=0,10000 output=test_atcorr rescale=0,10000 parameters=~/p6s.txt --v --o
(...)
Atmospheric correction...
 100%
Atmospheric correction complete.

r.info -r test_atcorr
min=0.5415085
max=8514.146

I get on Fedora 27 slightly different results:

r.info -r test_atcorr
min=0.5536021
max=8514.16

The 6s code is still numerically unstable. The 27,159 lines of code in i.atcorr need to be reviewed ...

This small patch makes sense and produces some visually reasonable output, therefore it could be applied to trunk and relbr74.

[mmetz]

Replying to mmetz:

This small patch makes sense and produces some visually reasonable output, therefore it could be applied to trunk and relbr74.

Applied to trunk and relbr74 in r72054,5.

[veroandreo]

Maybe relevant:

Updated spectral responses for S2: ​https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/document-library/-/asset_publisher/Wk0TKajiISaR/content/sentinel-2a-spectral-responses

[neteler]

Ticket retargeted after milestone closed

[mmetz]

Replying to mmetz:

Replying to mmetz:

This small patch makes sense and produces some visually reasonable output, therefore it could be applied to trunk and relbr74.

Applied to trunk and relbr74 in r72054,5.

Closing as fixed?

[sbl]

Yes, thanks!