Add AEC compression

.github/workflows/developer.yml

@@ -25,7 +25,7 @@ jobs:
     - name: install-dependencies
       run: |
         sudo apt-get update
-        sudo apt-get install libpng-dev zlib1g-dev libjpeg-dev doxygen gcovr valgrind
+        sudo apt-get install libpng-dev zlib1g-dev libjpeg-dev libaec-dev doxygen gcovr valgrind
 
     - name: "Build dependencies"
       uses: NOAA-EMC/ci-build-nceplibs@develop
@@ -36,7 +36,7 @@ jobs:
         w3emc-cmake-args: -DBUILD_WITH_BUFR=OFF
         ip-version: develop
         g2c-version: develop
-        g2c-cmake-args: -DBUILD_G2C=ON
+        g2c-cmake-args: -DBUILD_G2C=ON -DUSE_AEC=ON
 
     - name: checkout
       uses: actions/checkout@v4
@@ -57,7 +57,7 @@ jobs:
         mkdir build
         doxygen --version
         cd build
-        cmake -DBUILD_UTILS=OFF -DG2C_COMPARE=ON -DFTP_LARGE_TEST_FILES=ON -DENABLE_DOCS=ON -DCMAKE_PREFIX_PATH="$GITHUB_WORKSPACE/nceplibs/jasper;$GITHUB_WORKSPACE/nceplibs/NCEPLIBS-bacio;$GITHUB_WORKSPACE/nceplibs/NCEPLIBS-w3emc;$GITHUB_WORKSPACE/nceplibs/NCEPLIBS-ip;$GITHUB_WORKSPACE/nceplibs/NCEPLIBS-g2c" -DCMAKE_Fortran_FLAGS="-g -fprofile-abs-path -fprofile-arcs -ftest-coverage -O0 -Wall -fno-omit-frame-pointer -fsanitize=address" -DCMAKE_C_FLAGS="-g -fprofile-abs-path -fprofile-arcs -ftest-coverage -O0 -Wall -fno-omit-frame-pointer -fsanitize=address" -DFTP_TEST_FILES=ON -DTEST_FILE_DIR=/home/runner/data -DCMAKE_BUILD_TYPE=Debug ..
+        cmake -DBUILD_UTILS=OFF -DG2C_COMPARE=ON -DUSE_AEC=ON -DFTP_LARGE_TEST_FILES=ON -DENABLE_DOCS=ON -DCMAKE_PREFIX_PATH="$GITHUB_WORKSPACE/nceplibs/jasper;$GITHUB_WORKSPACE/nceplibs/NCEPLIBS-bacio;$GITHUB_WORKSPACE/nceplibs/NCEPLIBS-w3emc;$GITHUB_WORKSPACE/nceplibs/NCEPLIBS-ip;$GITHUB_WORKSPACE/nceplibs/NCEPLIBS-g2c" -DCMAKE_Fortran_FLAGS="-g -fprofile-abs-path -fprofile-arcs -ftest-coverage -O0 -Wall -fno-omit-frame-pointer -fsanitize=address" -DCMAKE_C_FLAGS="-g -fprofile-abs-path -fprofile-arcs -ftest-coverage -O0 -Wall -fno-omit-frame-pointer -fsanitize=address" -DFTP_TEST_FILES=ON -DTEST_FILE_DIR=/home/runner/data -DCMAKE_BUILD_TYPE=Debug ..
         make -j2 VERBOSE=1
 
     - name: test_asan

CMakeLists.txt

@@ -22,6 +22,7 @@ option(LOGGING "Turn on internal logging messages. Only useful to g2 developers.
 option(BUILD_4 "Build libg2_4.a" ON)
 option(BUILD_D "Build libg2_d.a" ON)
 option(BUILD_WITH_W3EMC "Build with NCEPLIBS-w3emc, enabling some GRIB1 functionality" ON)
+option(USE_AEC "Build with AEC (CCSDS) compression support" OFF)
 option(BUILD_UTILS "Build grib utilities" OFF)
 option(G2C_COMPARE "Enable copygb2 tests using g2c_compare" OFF)
 

spack/package.py

@@ -37,6 +37,7 @@ class G2(CMakePackage):
         when="@3.4.6:",
     )
     variant("w3emc", default=True, description="Enable GRIB1 through w3emc", when="@3.4.6:")
+    variant("aec", default=True, description="Use AEC library", when="@develop")
     variant("shared", default="False", when="@3.4.7:")
     variant("openmp", default=False, description="Use OpenMP multithreading")
     variant("utils", default=False, description="Build grib utilities")
@@ -45,6 +46,7 @@ class G2(CMakePackage):
     depends_on("jasper@:2.0.32", when="@:3.4.7")
     depends_on("jasper")
     depends_on("g2c", when="@develop")
+    depends_on("g2c@develop +aec", when="+aec")
     depends_on("libpng")
     depends_on("zlib-api")
     depends_on("bacio", when="@3.4.6:")
@@ -65,6 +67,7 @@ def cmake_args(self):
             self.define_from_variant("OPENMP", "openmp"),
             self.define_from_variant("CMAKE_POSITION_INDEPENDENT_CODE", "pic"),
             self.define_from_variant("BUILD_WITH_W3EMC", "w3emc"),
+            self.define_from_variant("USE_AEC", "aec"),
             self.define_from_variant("BUILD_SHARED_LIBS", "shared"),
             self.define("BUILD_4", self.spec.satisfies("precision=4")),
             self.define("BUILD_D", self.spec.satisfies("precision=d")),

src/CMakeLists.txt

@@ -15,6 +15,10 @@ if (BUILD_WITH_W3EMC)
   set(fortran_src ${fortran_src} gdt2gds.F90)
 endif()
 
+if (USE_AEC)
+  set(fortran_src ${fortran_src} g2aec.F90)
+endif()
+
 # These are the C source files.
 set(c_src mova2i.c)
 
@@ -42,6 +46,9 @@ foreach(kind ${kinds})
   # being built.
   target_compile_definitions(${lib_name}_f PUBLIC KIND=${kind})
 
+  if (USE_AEC)
+    target_compile_definitions(${lib_name}_f PRIVATE USE_AEC)
+  endif()
   # Build the C code.
   add_library(${lib_name}_c OBJECT ${c_src})
   target_include_directories(${lib_name}_c PRIVATE ${JASPER_INCLUDE_DIR})

src/g2aec.F90

@@ -0,0 +1,137 @@
+!> @file
+!> @brief Pack/unpack a data field that was packed with AEC compression.
+!> @author Eric Engle @date 2023-10-16
+
+!> Pack a data field into a AEC code stream as defined in
+!> [Data Representation Template
+!> 5.42](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-42.shtml).
+!>
+!> After the data are scaled, and the reference value is subtracted
+!> out, the data are passed to the AEC encoder.
+!> 
+!> This function also fills in GRIB2 Data Representation Template 5.42
+!> with the appropriate values.
+!>
+!> @param[in] fld The data values to pack.
+!> @param[in] width number of points in the x direction
+!> @param[in] height number of points in the y direction
+!> @param[inout] idrstmpl Contains the array of values for Data
+!> Representation Template [Table
+!> 5.42](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-42.shtml).
+!> - idrstmpl(1) Reference value - ignored on input.
+!> - idrstmpl(2) Binary Scale Factor.
+!> - idrstmpl(3) Decimal Scale Factor.
+!> - idrstmpl(4) Number of bits containing each grayscale pixel value
+!> - idrstmpl(5) Original field type, currently set = 0 on output
+!> Data values assumed to be reals.
+!> - idrstmpl(6) CCSDS compression options mask.
+!> - idrstmpl(7) Block size.
+!> - idrstmpl(8) Reference sample interval.
+!> @param[out] cpack The packed data field (character*1 array).
+!> @param[inout] lcpack When function is called, contains the length
+!> of buffer cpack. After functions returns, contains the length of
+!> the packed data in bytes.
+!>
+!> @author Eric Engle @date 2023-10-16
+subroutine aecpack(fld,width,height,idrstmpl,cpack,lcpack)
+
+  use, intrinsic :: iso_c_binding, only: c_size_t
+  implicit none
+
+  integer,intent(in) :: width,height
+  real,intent(in) :: fld(width*height)
+  character(len=1),intent(out) :: cpack(*)
+  integer,intent(inout) :: idrstmpl(*)
+  integer,intent(inout) :: lcpack
+  integer(c_size_t) :: width_c, height_c
+  integer :: ret
+
+  interface
+   function g2c_aecpackd(fld, width, height, idrstmpl, cpack, lcpack) bind(c)
+    use, intrinsic :: iso_c_binding
+    real(kind=c_double), intent(in):: fld(*)
+    integer(c_size_t), value, intent(in) :: width, height
+    integer(c_int), intent(inout) :: idrstmpl(*)
+    character(kind=c_char), intent(out) :: cpack(*)
+    integer(c_int), intent(out) :: lcpack
+    integer(c_int) :: g2c_aecpackd
+   end function g2c_aecpackd
+   function g2c_aecpackf(fld, width, height, idrstmpl, cpack, lcpack) bind(c)
+    use, intrinsic :: iso_c_binding
+    real(kind=c_float), intent(in):: fld(*)
+    integer(c_size_t), value, intent(in) :: width, height
+    integer(c_int), intent(inout) :: idrstmpl(*)
+    character(kind=c_char), intent(out) :: cpack(*)
+    integer(c_int), intent(out) :: lcpack
+    integer(c_int) :: g2c_aecpackf
+   end function g2c_aecpackf
+  end interface
+
+  width_c = width
+  height_c = height
+
+#if KIND==4
+  ret = g2c_aecpackf(fld, width_c, height_c, idrstmpl, cpack, lcpack)
+#else
+  ret = g2c_aecpackd(fld, width_c, height_c, idrstmpl, cpack, lcpack)
+#endif
+
+end subroutine
+
+!> Unpack a data field from a AEC code stream as defined in
+!> [Data Representation Template
+!> 5.42](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-42.shtml).
+!>
+!> @param[in] cpack The packed data field (character*1 array).
+!> @param[in] len length of packed field cpack().
+!> @param[in] idrstmpl Array of values for Data Representation
+!> [Template
+!> 5.42](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-42.shtml).
+!> @param[in] ndpts The number of data values to unpack.
+!> @param[out] fld Contains the unpacked data values.
+!>
+!> @author Eric Engle @date 2023-10-16
+subroutine aecunpack(cpack,len,idrstmpl,ndpts,fld)
+
+  use, intrinsic :: iso_c_binding, only: c_size_t
+  implicit none
+
+  character(len=1),intent(in) :: cpack(len)
+  integer,intent(in) :: ndpts,len
+  integer,intent(in) :: idrstmpl(*)
+  real,intent(out) :: fld(ndpts)
+  integer(c_size_t) :: ndpts_c, len_c
+  integer :: ret
+
+  interface
+   function g2c_aecunpackd(cpack, len, idrstmpl, ndpts, fld) bind(c)
+    use, intrinsic :: iso_c_binding
+    implicit none 
+    integer(c_size_t), value, intent(in) :: len
+    integer(c_size_t), value, intent(in) :: ndpts 
+    character(kind=c_char), intent(in) :: cpack(*)
+    integer(c_int), intent(in) :: idrstmpl(*)
+    real(kind=c_double), intent(out) :: fld(*)
+    integer(c_int) :: g2c_aecunpackd
+   end function g2c_aecunpackd
+   function g2c_aecunpackf(cpack, len, idrstmpl, ndpts, fld) bind(c)
+    use, intrinsic :: iso_c_binding
+    implicit none 
+    integer(c_size_t), value, intent(in) :: len
+    integer(c_size_t), value, intent(in) :: ndpts 
+    character(kind=c_char), intent(in) :: cpack(*)
+    integer(c_int), intent(in) :: idrstmpl(*)
+    real(kind=c_float), intent(out) :: fld(*)
+    integer(c_int) :: g2c_aecunpackf
+   end function g2c_aecunpackf
+  end interface
+
+  len_c = len
+  ndpts_c = ndpts
+#if KIND==4
+  ret = g2c_aecunpackf(cpack, len_c, idrstmpl, ndpts_c, fld)
+#else
+  ret = g2c_aecunpackd(cpack, len_c, idrstmpl, ndpts_c, fld)
+#endif
+
+end subroutine aecunpack

src/g2create.F90

@@ -509,6 +509,27 @@ end subroutine g2_sbytec1
      !print *, 'png size ', width, height
      call pngpack(pfld, width, height, idrstmpl, cpack, lcpack)
      !print *, 'png packed'
+#ifdef USE_AEC
+  elseif (idrsnum.eq.42) then ! AEC compression
+    if (ibmap.eq.255) then
+      call getdim(cgrib(lpos3), lensec3, width, height, iscan)
+      if (width.eq.0 .OR. height.eq.0) then
+          width=ndpts
+          height=1
+      elseif (width.eq.allones .OR. height.eq.allones) then
+          width=ndpts
+          height=1
+      elseif (ibits(iscan, 5, 1) .eq. 1) then ! Scanning mode: bit 3
+          itemp=width
+          width=height
+          height=itemp
+      endif
+    else
+      width=ndpts
+      height=1
+    endif
+    call aecpack(pfld, width, height, idrstmpl, cpack, lcpack);
+#endif /* USE_AEC */
   else
      print *, 'addfield: Data Representation Template 5.', idrsnum, &
           ' not yet implemented.'

src/g2get.F90

@@ -736,6 +736,12 @@ subroutine getfield(cgrib, lcgrib, ifldnum, igds, igdstmpl, &
            call pngunpack(cgrib(ipos + 5), lensec - 5, idrstmpl, &
                 ndpts, fld)
            have7 = .true.
+#ifdef USE_AEC
+        elseif (idrsnum .eq. 42) then
+           call aecunpack(cgrib(ipos + 5), lensec - 5, idrstmpl, &
+                ndpts, fld)
+            have7 = .true.
+#endif /* USE_AEC */
         else
            print *, 'getfield: Data Representation Template ', &
                 idrsnum, ' not yet implemented.'

src/g2jpc.F90

@@ -74,7 +74,6 @@ end function g2c_jpcpackf
   ret = g2c_jpcpackd(fld, width_c, height_c, idrstmpl, cpack, lcpack)
 #endif
 
-
 end subroutine
 
 !> Unpack a data field from a JPEG2000 code stream as defined in

src/g2unpack.F90

@@ -701,6 +701,10 @@ subroutine gf_unpack7(cgrib, lcgrib, iofst, igdsnum, igdstmpl,  &
      call jpcunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts, fld)
   elseif (idrsnum .eq. 41 .OR. idrsnum .eq. 40010) then
      call pngunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts, fld)
+#ifdef USE_AEC
+  elseif (idrsnum .eq. 42) then
+     call aecunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts, fld)
+#endif /* USE_AEC */
   else
      print *, 'gf_unpack7: Data Representation Template ', idrsnum, ' not yet implemented.'
      ierr = 4

tests/CMakeLists.txt

@@ -215,6 +215,10 @@ foreach(kind ${kinds})
     create_test(test_gdt2gds ${kind})
   endif()
 
+  if (USE_AEC)
+    create_test(test_aecpack ${kind})
+  endif()
+
   if (${kind} EQUAL "4")
     # Add this script to ensure that jasper is not issuing warnings.
     gu_test(run_jasper_warning_test)  

tests/test_aecpack.F90

@@ -0,0 +1,89 @@
+! This program tests the aecpack subroutine of the NCEPLIBS-g2
+! project.
+!
+! Alyson Stahl 2024-09-04
+program test_aecpack
+
+  implicit none
+
+  integer, parameter :: width=2, height=2, ndpts=4
+  real, parameter :: delta = 0.0000001
+  real :: fld(ndpts), fld2(ndpts)
+  integer :: idrstmpl(8)
+  character(len=1) :: cpack(80)
+  integer :: lcpack = 80
+  integer :: i
+
+  print *, 'Testing aecpack()/aecunpack()...'
+
+  ! Create the fld variable with data to pack
+  fld = (/1.0, 1.0, 1.0, 1.0/)
+  fld2 = (/0, 0, 0, 0/)
+
+  idrstmpl(1) = 0
+  idrstmpl(2) = 1
+  idrstmpl(3) = 1
+  idrstmpl(4) = 0
+  idrstmpl(5) = 0
+  idrstmpl(6) = 0
+  idrstmpl(7) = 16
+  idrstmpl(8) = 128
+
+  print *, 'Testing aecpack/aecunpack with no data...'
+
+  ! Pack the data.
+  call aecpack(fld, 0, height, idrstmpl, cpack, lcpack)
+  print *, 'lcpack: ', lcpack
+  if (lcpack .ne. 0) stop 1
+
+  ! Unpack the data.
+  call aecunpack(cpack, lcpack, idrstmpl, 0, fld2)
+
+  do i=1,ndpts
+    if (fld2(i) .ge. delta) then 
+      print *, fld2(i), ' does not equal zero'
+      stop 2
+    end if 
+  end do
+
+  print *, 'ok!'
+  print *, 'Testing aecpack/aecunpack with constant data...'
+
+  call aecpack(fld, width, height, idrstmpl, cpack, lcpack)
+  if (lcpack .ne. 0) stop 3
+
+  ! Unpack the data.
+  call aecunpack(cpack, lcpack, idrstmpl, ndpts, fld2)
+
+    ! Compare each value to see match, reals do not compare well
+  do i = 1, ndpts
+    if (abs(fld(i) - fld2(i)) .ge. delta) then
+       print *, fld(i), fld2(i), 'do not match'
+       stop 4
+    end if
+  end do
+
+  print *, 'ok!'
+  print *, 'Testing aecpack/aecunpack with data...'
+
+  fld = (/1.0, 2.0, 3.0, 4.0/)
+
+  ! Pack the data.
+  lcpack = 80
+  call aecpack(fld, width, height, idrstmpl, cpack, lcpack)
+  if (lcpack .le. 0) stop 5
+
+    ! Unpack the data.
+  call aecunpack(cpack, lcpack, idrstmpl, ndpts, fld2)
+
+  ! Compare each value to see match, reals do not compare well
+  do i = 1, ndpts
+     if (abs(fld(i) - fld2(i)) .ge. delta) then
+        print *, fld(i), fld2(i), 'do not match'
+        stop 6
+     end if
+  end do
+
+  print *, 'SUCCESS!'
+
+end program test_aecpack