Problem opening JP2 image using GDAL (interaction of -stats and -hist flags)

Question

I came across unexpected GDAL behaviour while working with Sentinel-2 data granules. The Sentinel data should have values between 0 and 10000 when opened correctly.

However, I have noticed that when working with the data, pixel values above 10000 are retrieved. After doing some digging, I found this post which points to similar problems in ESA's Snap software. They refer to a bug in the OpenJpeg library.

This could be the problem in my case as well, as my GDAL binaries use JP2OpenJPEG in the background (see below).

The strange behavior is the following: The values are in the expected range when using the gdalinfo -hist command, but are wrong when using gdalinfo -hist -stats.

I also get values above 10000 when using the GDAL c-api directly.

My questions are

1. Why is this happening?
2. Is there a way to open such jp2 images with 15bpp correctly in GDAL?

Here is an example and a link to the corresponding file:

http://sentinel-s2-l1c.s3.amazonaws.com/tiles/48/M/WA/2016/6/29/0/B11.jp2

$ gdalinfo B11.jp2 -hist
Band 1 Block=640x640 Type=UInt16, ColorInterp=Gray
0...10...20...30...40...50...60...70...80...90...100 - done.
  256 buckets from -15.7549 to 8050.75:
  23286152 2037325 ... 11251 
  Overviews: 2745x2745, 1372x1372, 686x686, 343x343
  Overviews: arbitrary
  Metadata:
    STATISTICS_MAXIMUM=8035
    STATISTICS_MEAN=274.9824
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=944.69220833573
  Image Structure Metadata:
    COMPRESSION=JPEG2000
    NBITS=15

$ rm B11.jp2.aux.xml
$ gdalinfo B11.jp2 -hist -stats
Band 1 Block=640x640 Type=UInt16, ColorInterp=Gray
  Minimum=0.000, Maximum=12356.000, Mean=272.236, StdDev=901.137
0...10...20...30...40...50...60...70...80...90...100 - done.
  256 buckets from -24.2275 to 12380.2:
  23286928 2372470 ... 0 0 1 0 1 
  Overviews: 2745x2745, 1372x1372, 686x686, 343x343
  Overviews: arbitrary
  Metadata:
    STATISTICS_MAXIMUM=12356
    STATISTICS_MEAN=272.236306847
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=901.13698170459
  Image Structure Metadata:
    COMPRESSION=JPEG2000
    NBITS=15

I am using Ubuntu 16.06 with the Ubuntugis gdal version

$ gdalinfo --version
GDAL 2.1.0, released 2016/04/25

$ gdalinfo --formats | grep JP2
JP2OpenJPEG -raster,vector- (rwv): JPEG-2000 driver based on OpenJPEG library

Answer 1

A test with GDAL 2.3-dev and OpenJPEG driver gives different result with -stats:

Band 1 Block=640x640 Type=UInt16, ColorInterp=Gray
  Minimum=0.000, Maximum=12356.000, Mean=272.236, StdDev=901.137
  Overviews: 2745x2745, 1372x1372, 686x686, 343x343
  Overviews: arbitrary
  Metadata:
    STATISTICS_MAXIMUM=12356
    STATISTICS_MEAN=272.236306847
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=901.13698170459
  Image Structure Metadata:
    COMPRESSION=JPEG2000
    NBITS=15

I am getting also the same result with Kakadu driver:

Driver: JP2KAK/JPEG-2000 (based on Kakadu v7.3.2)
Files: b11.jp2
...
Metadata:
  STATISTICS_MAXIMUM=12356
  STATISTICS_MEAN=272.236306847
  STATISTICS_MINIMUM=0
  STATISTICS_STDDEV=901.13698170459

Are you sure that the file does not contain values higher than 10000? With which method have you confirmed it?

Answer 2

I found out that the band values above 10'000 make sense after all.

10k is the number that represents 100% reflectance in the data. That's why I thought that the values should be bound by this number, but they are not.

Reflectances can be higher than 100%, see the answer of a related question on the step forum.

I think the difference in the output of GDAL might be related to approximate vs excact calculations of the statistics.