GDAL on Windows/Mac do not properly decode GRIB2 compressed with the JPEG2000

archambaultb commented 3 years ago

Expected behavior and actual behavior.

We discovered that on Windows and Mac, some versions of GDAL do not properly decode GRIB2 files if the data is compressed with the JPEG2000 option. When GDAL decodes the data, the values are orders of magnitude off (MAX values around 10^5 instead of 10^2 for the APCP parameter, precipitation analysis). The problem does not happen on the Linux platform neither when using other tools to read the data (such as wgrib2, on Linux).

Actual behavior: GRIB2 compressed with the JPEG2000 decoded on Windows/Mac gives values orders of magnitude off.
Expected behavior: GRIB2 compressed with the JPEG2000 decoded properly, like it's done under correctly under Linux.

Steps to reproduce the problem.

First let’s get the real min/max values of the GRIB2 sample file on a Linux platform with gdalinfo:

Version tested:

Linux: Ubuntu 18.04.2 LTS (Bionic Beaver)
GDAL/OGR GDAL 2.3.2, released 2018/09/21
Note all versions tested of Linux/GDAL gives the same result.

& gdalinfo -mm https://collaboration.cmc.ec.gc.ca/cmc/cmos/gdal_Jpeg2000/CMC_RDPA_APCP-024-0700cutoff_SFC_0_ps10km_2021041106_000.grib2 | grep "Min/Max"
    Computed Min/Max=0.000,122.170
    Computed Min/Max=0.000,0.916

This is the expected result:

122.170 is the max value of the precipitation in kg/m^2.
0.916 is the max value of the confidence index (defined between 0 and 1)

Now if we try to get the min/max values of the same GRIB2 sample file on a Windows platform with gdalinfo:

Version tested:

Windows: Windows 10 Enterprise
QGIS: 3.10.2-A Coruña
GDAL/OGR 3.0.3 (installed within QGIS)

Steps to reproduce the bug.

Open Windows Powershell
Move to the bin directory of your QGIS installation

cd 'C:\Program Files (x86)\QGIS 3.10\bin'

Get Max/min of the GRIB2 sample file using gdalinfo

./gdalinfo.exe -mm https://collaboration.cmc.ec.gc.ca/cmc/cmos/gdal_Jpeg2000/CMC_RDPA_APCP-024-0700cutoff_SFC_0_ps10km_2021041106_000.grib2 | Select-String -Pattern "Min/Max"
Computed Min/Max=0.000,**268435.469**
Computed Min/Max=0.000,**2147.484**

268435.469 is definitely not a realistic value for a Precipitation Max in kg/m^2

And 2147.484 is not possible for a parameter defined between 0 and 1.

Get Max/min of the GRIB2 sample file repacked using another compression than JPEG200 using gdalinfo (simple packing is used for the example)

./gdalinfo.exe -mm https://collaboration.cmc.ec.gc.ca/cmc/cmos/gdal_Jpeg2000/CMC_RDPA_APCP-024-0700cutoff_SFC_0_ps10km_2021041106_000_simple_packing.grib2 | Select-String -Pattern "Min/Max"
Computed Min/Max=0.000,122.170
Computed Min/Max=0.000,0.916

Which gives the exact same value as the gdalinfo calls on Linux.

Since the only difference between the two sample files is the compression algorithm, it seems that there is a bug in the JPEG2000 decoding on GDAL 3.0.3 (Windows).

Other versions tested

Here is some other versions and configuration that we have tested:

OS	GIS	GDAL	Bug?
Windows	QGIS 3.16.3-Hannover	3.1.4	Bug is present
Windows 10	QGIS 3.10.2-A Coruna	3.0.3	Bug is present
Windows 10	ArcGIS (Explorer desktop)	3.0.3	Bug is present
Windows 10	QGIS 3.10.1-A Coruna	3.0.2	Bug is present
Windows	QGIS 2.18.28	2.4.0	Bug is present
Mac	QGIS 3.18.0	3.2.1	Bug is present
Linux Ubuntu 18.04	QGIS 3.18.2	3.2.2	No problems
Linux Ubuntu 18.04	QGIS 2	3.0.3	No problems
Windows 10	No GIS used, GDAL tested within Conda	2.4.3 2.4.4 3.1.4 3.2.0 3.2.1 3.2.2	No problems
Windows 10	No GIS used, GDAL Internals https://www.gisinternals.com/release.php	3.2.0	No problems

As you can see, we have tested the same sample data with GDAL on Windows via Conda and there was no issue. This led us to believe it may be related to the way QGIS and ArcGIS integrates the GDAL GRIB2 driver or the JPEG2000 compression library."

Operating system

See the bug description and the table above.

GDAL version and provenance

See the bug description and the table above.

rouault commented 3 years ago

Please edit your above table by mentioning with JPEG2000 is actually used (this is the first line of the gdalinfo output, like JP2OpenJPEG, JP2ECW, etc.), and which version of the corresponding library (openjpeg, ECW SDK, etc.) is used.

archambaultb commented 3 years ago

Hi, thank @rouault you for the quick feedback.

I have ran the gdalinfo on the first Windows case described in the ticket (QGIS: 3.10.2-A Coruña - GDAL/OGR 3.0.3 (installed within QGIS)).

First line indicates Driver: GRIB/GRIdded Binary (.grb, .grb2) is this the info you need? Before updating the table, I want to make sure that I have the right info.

Here is the full gdalinfo output I have for QGIS: 3.10.2-A Coruña - GDAL/OGR 3.0.3 (installed within QGIS) :

PS C:\Program Files (x86)\QGIS 3.10\bin> .\gdalinfo.exe -mm https://collaboration.cmc.ec.gc.ca/cmc/cmos/gdal_Jpeg2000/CMC_RDPA_APCP-024-0700cutoff_SFC_0_ps10km_2021041106_000.grib2
Driver: GRIB/GRIdded Binary (.grb, .grb2)
Files: none associated
Size is 935, 824
Coordinate System is:
PROJCRS["unnamed",
    BASEGEOGCRS["Coordinate System imported from GRIB file",
        DATUM["unnamed",
            ELLIPSOID["Sphere",6371229,0,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["unnamed",
        METHOD["Polar Stereographic (variant B)",
            ID["EPSG",9829]],
        PARAMETER["Latitude of standard parallel",60,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8832]],
        PARAMETER["Longitude of origin",249,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8833]],
        PARAMETER["False easting",0,
            LENGTHUNIT["Metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["Metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",south,
            ORDER[1],
            LENGTHUNIT["Metre",1]],
        AXIS["northing",south,
            ORDER[2],
            LENGTHUNIT["Metre",1]]]
Data axis to CRS axis mapping: 1,2
Origin = (-4556441.403315245173872,920682.141165950335562)
Pixel Size = (10000.000000000000000,-10000.000000000000000)
Corner Coordinates:
Upper Left  (-4556441.403,  920682.141) (147d34'35.51"E, 47d17'21.08"N)
Lower Left  (-4556441.403,-7319317.859) (142d54'11.52"W, 18d 6' 5.46"N)
Upper Right ( 4793558.597,  920682.141) ( 10d 7'40.08"W, 45d21'25.84"N)
Lower Right ( 4793558.597,-7319317.859) ( 77d46'42.45"W, 17d17'57.97"N)
Center      (  118558.597,-3199317.859) (108d52'39.85"W, 59d51'25.44"N)
Band 1 Block=935x1 Type=Float64, ColorInterp=Undefined
Band 1 Block=935x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] SFC="Ground or water surface"
    Computed Min/Max=0.000,268435.469
  NoData Value=9999
  Metadata:
    GRIB_COMMENT=Total precipitation [kg/(m^2)]
    GRIB_DISCIPLINE=0(Meteorological)
    GRIB_ELEMENT=APCP
    GRIB_FORECAST_SECONDS=86400 sec
    GRIB_IDS=CENTER=54(Montreal) SUBCENTER=0 MASTER_TABLE=4 LOCAL_TABLE=0 SIGNF_REF_TIME=0(Analysis) REF_TIME=2021-04-11T06:00:00Z PROD_STATUS=0(Operational) TYPE=0(Analysis)
    GRIB_PDS_PDTN=8
    GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES=1 8 0 30 30 0 0 1 24 1 0 0 255 -127 -2147483647 2021 4 11 6 0 0 1 0 1 2 1 4294967272 1 0
    GRIB_PDS_TEMPLATE_NUMBERS=1 8 0 30 30 0 0 0 1 0 0 0 24 1 0 0 0 0 0 255 255 255 255 255 255 7 229 4 11 6 0 0 1 0 0 0 0 1 2 1 255 255 255 232 1 0 0 0 0
    GRIB_REF_TIME=  1618120800 sec UTC
    GRIB_SHORT_NAME=0-SFC
    GRIB_UNIT=[kg/(m^2)]
    GRIB_VALID_TIME=  1618120800 sec UTC
Band 2 Block=935x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] SFC="Ground or water surface"
    Computed Min/Max=0.000,2147.484
  NoData Value=9999
  Metadata:
    GRIB_COMMENT=(prodType 0, cat 1, subcat 193) [-]
    GRIB_DISCIPLINE=0(Meteorological)
    GRIB_ELEMENT=unknown
    GRIB_FORECAST_SECONDS=86400 sec
    GRIB_IDS=CENTER=54(Montreal) SUBCENTER=0 MASTER_TABLE=4 LOCAL_TABLE=1 SIGNF_REF_TIME=0(Analysis) REF_TIME=2021-04-11T06:00:00Z PROD_STATUS=0(Operational) TYPE=0(Analysis)
    GRIB_PDS_PDTN=8
    GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES=1 193 0 30 30 0 0 1 24 1 0 0 255 -127 -2147483647 2021 4 11 6 0 0 1 0 0 2 1 4294967272 1 0
    GRIB_PDS_TEMPLATE_NUMBERS=1 193 0 30 30 0 0 0 1 0 0 0 24 1 0 0 0 0 0 255 255 255 255 255 255 7 229 4 11 6 0 0 1 0 0 0 0 0 2 1 255 255 255 232 1 0 0 0 0
    GRIB_REF_TIME=  1618120800 sec UTC
    GRIB_SHORT_NAME=0-SFC
    GRIB_UNIT=[-]
    GRIB_VALID_TIME=  1618120800 sec UTC

rouault commented 3 years ago

First line indicates Driver: GRIB/GRIdded Binary (.grb, .grb2) is this the info you need?

ah sorry, I didn't think that it was grib wrapping jpeg2000, and not plain jpeg2000. The easiest is probably that you do "gdalinfo some_random_jpeg2000_file.jp2" where such random file can be https://github.com/OSGeo/gdal/blob/master/autotest/gdrivers/data/jpeg2000/byte.jp2?raw=true to see which active jpeg2000 driver you have And you can confirm that the JPEG2000 driver reported with that JP2 file is actually used by the GRIB driver by doing "gdalinfo --debug on -mm your.grb" and checking a line like "GDAL: GDALOpen(/vsimem/workgrib{some_hexadecimal_number}.jpc, this={some_hexadecimal_number}) succeeds as {JPEG2000_DRIVER_NAME_HERE}"

archambaultb commented 3 years ago

I think the JPEG Driver used is JP2ECW:

Here is the info I got with your file jpeg2000/byte.jp2

PS C:\Program Files (x86)\QGIS 3.10\bin> .\gdalinfo.exe -mm  https://github.com/OSGeo/gdal/blob/master/autotest/gdrivers/data/jpeg2000/byte.jp2?raw=true
Driver: JP2ECW/ERDAS JPEG2000 (SDK 5.3)

And the gdalinfo --debug on my file gives this:

PS C:\Program Files (x86)\QGIS 3.10\bin> .\gdalinfo.exe --debug on -mm https://collaboration.cmc.ec.gc.ca/cmc/cmos/gdal_Jpeg2000/CMC_RDPA_APCP-024-0700cutoff_SFC_0_ps10km_2021041106_000.grib2
...
GDAL: GDALOpen(/vsimem/work_grib_0420D913.jpc, this=00C95848) succeeds as JP2ECW

rouault commented 3 years ago

@jef-n Which ECW SDK version is used with latest QGIS 3.18 ? in OSGeo4W

rouault commented 3 years ago

I can also reproduce the issue with ECW SDK 5.5 on Linux.

jef-n commented 3 years ago

@jef-n Which ECW SDK version is used with latest QGIS 3.18 ? in OSGeo4W

5.5 on testing 5.3 on otherwise

VirtualTim commented 1 year ago

The ECWSDK only supports up to 28 bit JP2 decoding/encoding. From the user guide (page 30)

The unified compression architecture of JPEG 2000 enables seamless image component compression from 1 to 28 bits deep.

rouault commented 1 year ago

The ECWSDK only supports up to 28 bit JP2 decoding/encoding.

That shouldn't be an issue. This is just 20 bit JP2 :

$ gdalinfo /vsisubfile/96499_,CMC_RDPA_APCP-024-0700cutoff_SFC_0_ps10km_2021041106_000.grib2
Driver: JP2OpenJPEG/JPEG-2000 driver based on OpenJPEG library
Files: /vsisubfile/96499_,CMC_RDPA_APCP-024-0700cutoff_SFC_0_ps10km_2021041106_000.grib2
Size is 750157, 1
Image Structure Metadata:
  COMPRESSION_REVERSIBILITY=LOSSLESS (possibly)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,    1.0)
Upper Right (750157.0,    0.0)
Lower Right (750157.0,    1.0)
Center      (375078.5,    0.5)
Band 1 Block=750157x1 Type=UInt32, ColorInterp=Undefined
  Overviews: 375079x1, 187540x1, 93770x1, 46885x1, 23443x1
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000
    NBITS=20

Perhaps the issue is that it is a big one-line JPEG2000 ? (GRIB2 uses such a weird formulation)

OSGeo / gdal