GRIB: implement CCSDS Adaptive Entropy Coding decompression (template…

… 5.42).

Requires building GDAL against libaec (Cf https://gitlab.dkrz.de/k202009/libaec)

Cf NCAR/ncl#204

Fixes #8092

.github/workflows/cmake_builds.yml

@@ -54,7 +54,7 @@ jobs:
       run: |
         sudo add-apt-repository -y ppa:ubuntugis/ubuntugis-unstable
         sudo apt-get update
-        sudo apt-get install -y -q bison libjpeg-dev libgif-dev liblzma-dev libzstd-dev libgeos-dev git \
+        sudo apt-get install -y -q bison libaec-dev libjpeg-dev libgif-dev liblzma-dev libzstd-dev libgeos-dev git \
            libcurl4-gnutls-dev libproj-dev libxml2-dev  libxerces-c-dev libnetcdf-dev netcdf-bin \
            libpoppler-dev libpoppler-private-dev gpsbabel libhdf4-alt-dev libhdf5-serial-dev libpodofo-dev poppler-utils \
            libfreexl-dev unixodbc-dev libwebp-dev libepsilon-dev liblcms2-2 libcrypto++-dev libkml-dev \
@@ -430,7 +430,7 @@ jobs:
       shell: bash -l {0}
       run: |
         conda install --yes --quiet curl libiconv icu python=3.10 swig numpy pytest pytest-env filelock zlib lxml jsonschema
-        conda install --yes --quiet proj geos hdf4 hdf5 kealib \
+        conda install --yes --quiet libaec proj geos hdf4 hdf5 kealib \
             libnetcdf openjpeg poppler libtiff libpng xerces-c expat libxml2 kealib json-c \
             cfitsio freexl geotiff libjpeg-turbo libpq libspatialite libwebp-base pcre pcre2 postgresql \
             sqlite tiledb zstd cryptopp cgal doxygen librttopo libkml openssl xz \

.github/workflows/ubuntu_20.04/Dockerfile.ci

@@ -21,6 +21,7 @@ RUN apt-get update -y \
     git \
     gnupg \
     lcov \
+    libaec-dev \
     libarmadillo-dev \
     libblosc-dev \
     libboost-dev \

autotest/gdrivers/grib.py

@@ -2230,3 +2230,15 @@ def test_grib_grib2_wrong_earth_shape():
         ds = gdal.Open("data/grib/byte_wrong_earth_shape.grib2")
         assert ds.GetSpatialRef().GetSemiMajor() == 6377563.396
         assert ds.GetSpatialRef().GetSemiMinor() == 6377563.396
+
+
+# Test reading file with template 5.42 / CCSDS szip/aes compression
+
+
+def test_grib_grib2_template_5_42_CCDS_aes_decompression():
+
+    ds = gdal.Open("data/grib/template_5_42_ccsds_aec.grb2")
+    if gdal.GetDriverByName("GRIB").GetMetadataItem("HAVE_AEC"):
+        assert ds.GetRasterBand(1).Checksum() == 41970
+    else:
+        assert ds.GetRasterBand(1).Checksum() == -1

cmake/helpers/CheckDependentLibraries.cmake

@@ -713,6 +713,9 @@ gdal_check_package(Blosc "Blosc compression" CAN_DISABLE)
 define_find_package2(ARCHIVE archive.h archive)
 gdal_check_package(ARCHIVE "Multi-format archive and compression library library (used for /vsi7z/" CAN_DISABLE)
 
+define_find_package2(LIBAEC libaec.h aec)
+gdal_check_package(LIBAEC "Adaptive Entropy Coding implementing Golomb-Rice algorithm (used by GRIB)" CAN_DISABLE)
+
 define_find_package2(JXL jxl/decode.h jxl PKGCONFIG_NAME libjxl)
 gdal_check_package(JXL "JPEG-XL compression" CAN_DISABLE)
 

doc/source/development/building_from_source.rst

@@ -1026,6 +1026,24 @@ of the original input image is preserved (within user defined error bounds).
     Control whether to use the LERC internal library. Defaults depends on GDAL_USE_INTERNAL_LIBS. When set
     to ON, has precedence over GDAL_USE_LERC=ON
 
+LIBAEC
+******
+
+`libaec <https://gitlab.dkrz.de/k202009/libaec>`_ is a compression library which offers
+the extended Golomb-Rice coding as defined in the CCSDS recommended standard 121.0-B-3.
+It is used by the :ref:`raster.grib` driver.
+
+.. option:: LIBAEC_INCLUDE_DIR
+
+    Path to an include directory with the ``libaec.h`` header file.
+
+.. option:: LIBAEC_LIBRARY
+
+    Path to a shared or static library file.
+
+.. option:: GDAL_USE_LIBAEC=ON/OFF
+
+    Control whether to use LIBAEC. Defaults to ON when LIBAEC is found.
 
 LibKML
 ******

docker/alpine-normal/Dockerfile

@@ -44,24 +44,9 @@ RUN apk add --no-cache \
     apache-arrow-dev \
     # For spatialite (and GDAL)
     libxml2-dev \
+    libaec-dev \
     && mkdir -p /build_thirdparty/usr/lib
 
-# Build szip
-ARG SZIP_VERSION=2.1.1
-RUN if test "${SZIP_VERSION}" != ""; then ( \
-    wget -q https://support.hdfgroup.org/ftp/lib-external/szip/${SZIP_VERSION}/src/szip-${SZIP_VERSION}.tar.gz \
-    && tar xzf szip-${SZIP_VERSION}.tar.gz \
-    && rm -f szip-${SZIP_VERSION}.tar.gz \
-    && cd szip-${SZIP_VERSION} \
-    && CFLAGS=-O2 ./configure --prefix=/usr --disable-static \
-    && make -j$(nproc) \
-    && make install \
-    && cp -P /usr/lib/libsz*.so* /build_thirdparty/usr/lib \
-    && for i in /build_thirdparty/usr/lib/*; do strip -s $i 2>/dev/null || /bin/true; done \
-    && cd .. \
-    && rm -rf szip-${SZIP_VERSION} \
-    ); fi
-
 # Build hdf4
 ARG HDF4_VERSION=4.2.16
 RUN if test "${HDF4_VERSION}" != "" -a "$(uname -m)" = "x86_64"; then ( \
@@ -285,6 +270,7 @@ RUN apk add --no-cache \
         unixodbc \
         libpq \
         apache-arrow \
+        libaec \
     # Remove /usr/lib/libopenjp2.so.2.3.0 since we are building v2.3.1 manually
     # && rm -f /usr/lib/libopenjp2.so.2.3.0 \
     # libturbojpeg.so is not used by GDAL. Only libjpeg.so*

docker/ubuntu-full/Dockerfile

@@ -67,6 +67,7 @@ RUN . /buildscripts/bh-set-envvars.sh \
        libdeflate-dev${APT_ARCH_SUFFIX} libblosc-dev${APT_ARCH_SUFFIX} liblz4-dev${APT_ARCH_SUFFIX} libbz2-dev${APT_ARCH_SUFFIX} \
        libbrotli-dev${APT_ARCH_SUFFIX} \
        libarchive-dev${APT_ARCH_SUFFIX} \
+       libaec-dev${APT_ARCH_SUFFIX}
     && rm -rf /var/lib/apt/lists/*
 
 # Build likbkea
@@ -310,6 +311,7 @@ RUN apt-get update \
         libdeflate0 libblosc1 liblz4-1 \
         libbrotli1 \
         libarchive13 \
+        libaec0 \
         python-is-python3 \
         # pil for antialias option of gdal2tiles
         python3-pil \

frmts/grib/CMakeLists.txt

@@ -122,3 +122,9 @@ endif ()
 if (GDAL_USE_ZLIB AND NOT GDAL_USE_ZLIB_INTERNAL)
   gdal_target_link_libraries(gdal_GRIB PRIVATE ZLIB::ZLIB)
 endif ()
+
+if (GDAL_USE_LIBAEC)
+  target_compile_definitions(gdal_GRIB PRIVATE -DUSE_AEC)
+  gdal_target_link_libraries(gdal_GRIB PRIVATE LIBAEC::LIBAEC)
+  target_sources(gdal_GRIB PRIVATE degrib/g2clib/aecunpack.c)
+endif ()

frmts/grib/degrib/g2clib/aecunpack.c

@@ -0,0 +1,82 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <limits.h>
+#include "grib2.h"
+
+#include "libaec.h"
+
+// Cf https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-42.shtml
+// and https://github.com/erget/wgrib2/commit/07b0f6fcb9669e0e3285318f50d516731b6956b2
+
+g2int aecunpack(unsigned char *cpack,g2int len,g2int *idrstmpl,g2int ndpts,
+                g2float *fld)
+{
+
+      g2int  iret = 0;
+      g2float  refV;
+
+      rdieee(idrstmpl+0,&refV,1);
+      g2float bscale = (g2float)int_power(2.0,idrstmpl[1]);
+      g2float dscale = (g2float)int_power(10.0,-idrstmpl[2]);
+      g2float bdscale = bscale * dscale;
+      g2float refD = refV * dscale;
+
+      g2int nbits = idrstmpl[3];
+//
+//  if nbits equals 0, we have a constant field where the reference value
+//  is the data value at each gridpoint
+//
+      if (nbits != 0) {
+         int nbytes_per_sample = (nbits + 7) / 8;
+         if( ndpts != 0 && nbytes_per_sample > INT_MAX / ndpts )
+         {
+             return 1;
+         }
+         g2int* ifld=(g2int *)calloc(ndpts,sizeof(g2int));
+         // Was checked just before
+         // coverity[integer_overflow,overflow_sink]
+         unsigned char* ctemp=(unsigned char *)calloc(ndpts * nbytes_per_sample,1);
+         if ( ifld == NULL || ctemp == NULL) {
+            fprintf(stderr, "Could not allocate space in aecunpack.\n"
+                    "Data field NOT unpacked.\n");
+            free(ifld);
+            free(ctemp);
+            return(1);
+         }
+
+         struct aec_stream strm = {0};
+         strm.flags = idrstmpl[5]; // CCSDS compression options mask
+         strm.bits_per_sample = nbits;
+         strm.block_size = idrstmpl[6];
+         strm.rsi = idrstmpl[7]; // Restart interval
+
+         strm.next_in = cpack;
+         strm.avail_in = len;
+         strm.next_out = ctemp;
+         strm.avail_out = ndpts * nbytes_per_sample;
+
+         // Note: libaec doesn't seem to be very robust to invalid inputs...
+         int status = aec_buffer_decode(&strm);
+         if (status != AEC_OK)
+         {
+             fprintf(stderr, "aec_buffer_decode() failed with return code %d", status);
+             iret = 1;
+         }
+         else
+         {
+             gbits(ctemp,ndpts * nbytes_per_sample,ifld,0,nbytes_per_sample*8,0,ndpts);
+             g2int j;
+             for (j=0;j<ndpts;j++) {
+                fld[j] = refD + bdscale*(g2float)(ifld[j]);
+             }
+         }
+         free(ctemp);
+         free(ifld);
+      }
+      else {
+         g2int j;
+         for (j=0;j<ndpts;j++) fld[j]=refD;
+      }
+
+      return(iret);
+}

frmts/grib/degrib/g2clib/drstemplates.c

@@ -14,14 +14,16 @@ static const struct drstemplate templatesdrs[MAXDRSTEMP] = {
          { 4, 1, 0, {1} },
              // 5.50: Spectral Data - Simple Packing
          { 50, 5, 0, {4,-2,-2,1,4} },
-             // 5.51: Spherical Harmonics data - Complex packing 
+             // 5.51: Spherical Harmonics data - Complex packing
          { 51, 10, 0, {4,-2,-2,1,-4,2,2,2,4,1} },
 //           // 5.1: Matrix values at gridpoint - Simple packing
 //         { 1, 15, 1, {4,-2,-2,1,1,1,4,2,2,1,1,1,1,1,1} },
              // 5.40: Grid point data - JPEG2000 encoding
          { 40, 7, 0, {4,-2,-2,1,1,1,1} },
              // 5.41: Grid point data - PNG encoding
          { 41, 5, 0, {4,-2,-2,1,1} },
+             // 5.42: Grid point and spectral data - CCSDS szip
+         { 42, 8, 0, {4,-2,-2,1,1,1,1,-2} },
              // 5.40000: Grid point data - JPEG2000 encoding
          { 40000, 7, 0, {4,-2,-2,1,1,1,1} },
              // 5.40010: Grid point data - PNG encoding

frmts/grib/degrib/g2clib/drstemplates.h

@@ -31,7 +31,7 @@
 //
 ///////////////////////////////////////////////////////////////////////
 
-      #define MAXDRSTEMP 10             // maximum number of templates
+      #define MAXDRSTEMP 11             // maximum number of templates
       #define MAXDRSMAPLEN 200          // maximum template map length
 
       struct drstemplate

frmts/grib/degrib/g2clib/g2_unpack7.c

@@ -5,10 +5,6 @@
 #include <string.h>
 #include "grib2.h"
 
-#ifdef USE_PNG
-  g2int pngunpack(unsigned char *,g2int,g2int *,g2int, g2float *);
-#endif  /* USE_PNG */
-
 static float DoubleToFloatClamp(double val) {
    if (val >= FLT_MAX) return FLT_MAX;
    if (val <= -FLT_MAX) return -FLT_MAX;
@@ -216,6 +212,17 @@ g2int g2_unpack7(unsigned char *cgrib,g2int cgrib_length,g2int *iofst,g2int igds
         pngunpack(cgrib+ipos,lensec-5,idrstmpl,ndpts,lfld);
         }
 #endif  /* USE_PNG */
+      else if (idrsnum == 42) {
+#ifdef USE_AEC
+        aecunpack(cgrib+ipos,lensec-5,idrstmpl,ndpts,lfld);
+#else
+        fprintf(stderr,"g2_unpack7: Data Representation Template 5.42 decoding requires building against libaec.\n");
+        ierr=4;
+        if ( lfld != 0 ) free(lfld);
+        *fld=0;     //NULL
+        return(ierr);
+#endif
+        }
       else {
         fprintf(stderr,"g2_unpack7: Data Representation Template 5.%d not yet implemented.\n",(int)idrsnum);
         ierr=4;

frmts/grib/degrib/g2clib/gdal_g2clib_symbol_rename.h

@@ -1,4 +1,5 @@
 /* This is a generated file by rename_g2clib_symbols.h. *DO NOT EDIT MANUALLY !* */
+#define aecunpack gdal_aecunpack
 #define cmplxpack gdal_cmplxpack
 #define compack gdal_compack
 #define comunpack gdal_comunpack

frmts/grib/degrib/g2clib/grib2.h

@@ -14,7 +14,7 @@
 // 2009-01-14  Vuong     Changed struct template to gtemplate
 //
 //   Each element of structure gribfield is defined as:
-//   
+//
 //   gribfield gfld;
 //
 //        gfld->version = GRIB edition number ( currently 2 )
@@ -294,5 +294,8 @@ void pngpack(g2float *fld,g2int width,g2int height,g2int *idrstmpl,
 
 int dec_jpeg2000(const void *injpc,g2int bufsize,g2int **outfld,g2int outpixels);
 
+g2int aecunpack(unsigned char *cpack,g2int len,g2int *idrstmpl,g2int ndpts,
+                g2float *fld);
+
 #endif  /*  grib2_H  */
 

frmts/grib/gribdataset.cpp

@@ -2920,5 +2920,9 @@ void GDALRegister_GRIB()
     poDriver->pfnCreateCopy = GRIBDataset::CreateCopy;
     poDriver->pfnUnloadDriver = GDALDeregister_GRIB;
 
+#ifdef USE_AEC
+    poDriver->SetMetadataItem("HAVE_AEC", "YES");
+#endif
+
     GetGDALDriverManager()->RegisterDriver(poDriver);
 }