Merge pull request #298 from NOAA-EMC/c3_24bit_overflow_WNE

issue #292 fix overflow for complex-3 packing

docs/web_docs/set_grib_max_bits.html

@@ -133,6 +133,30 @@ <h3>Introduction</h3>
 to 25 after accounting for the limits of the packing and
 unpacking routines.
 
+<p> The complex-2 and complex-3 compression schemes also
+store scaled integers but different.  The complex-1 stores
+the scaled intger values, so it is limited to a 25 bit
+integers.  Hewever, complex-2 stores the delta (the difference
+from previously defined grid point).  The magnitude of delta
+is limited by the original series but the delta can be
+either positive or negative.  So under worst case, the
+deltas need one more bit to store the scaled integers.
+So if you scale the grid point values to 24 bits, you
+will have a scaled deltas being 25 bits.  The complex-3
+packing stores the deltas of the delta series. This is
+good for smooth fields, but for the worst case could 
+cause two extra bits to store the original grid data.
+So at worst case, the data would need to be scaled at 23
+bits.
+
+<p> The extra bits needed for complex-2 and complex-3 are
+very rare.  It was first encountered in 11/2024 which
+is about a decade after the complex packing was introduced.
+In 11/2024, code was added to the complex packing to detect
+the increase in bits needed and reduce the precision if it
+cause the precision to be greater than 25 bits.
+
+
 <p> Grib fields are usually stored as a scaled integers,
 and usually you don't need 25 bits of precision.  For
 example, temperature to the nearest 0.1 degree, means

tests/CMakeLists.txt

@@ -109,6 +109,7 @@ copy_test_data(ref_simple_packing.grib2.spread.txt)
 copy_test_data(ref_new_grid_gdt_32769.grib2)
 copy_test_data(png_4bits.png)
 copy_test_data(large_png.grb2)
+copy_test_data(ref_c3_overflow.grib2)
 
 # Run these shell tests.
 shell_test(run_wgrib2_tests)

tests/run_wgrib2_tests.sh

@@ -145,7 +145,7 @@ if [ `echo "$line" | grep -c ":rpn_rms=0:"` -ne 1 ] ; then
   exit 1
 fi
 
-
+../wgrib2/wgrib2 data/ref_c3_overflow.grib2 -set_grib_type c3 -set_grib_max_bits 25 -grib_out c3_overflow.grib2
 
 echo "*** SUCCESS!"
 exit 0

wgrib2/complex_pk.c

@@ -84,13 +84,13 @@ static int sizeofsection2(int mn, int mx, int n, int ref_bits, int width_bits, i
 }
 
 static int size_all(struct section *s, int ref_bits, int width_bits, int has_undef) {
-	unsigned int bytes, bits;
+        unsigned int bytes, bits;
 
         bytes = bits = 0;
         while (s) {
             bits += sizeofsection(s, ref_bits, width_bits, has_undef);
-	    bytes += (bits / 8);
-	    bits = bits % 8;
+            bytes += (bits / 8);
+            bits = bits % 8;
             s = s->tail;
         }
         return (int) (bytes + (bits != 0));
@@ -722,8 +722,7 @@ int complex_grib_out(unsigned char **sec, float *data, unsigned int ndata,
 	if (int_min_max_array(v, nndata, &vmn, &vmx) != 0) fatal_error("complex_pk: int_min_max error","");
     }
 
-//fprintf(stderr , "\n pre process v[i] data extri_0 %d extra_1 %d\n",extra_0, extra_1);
-//for (i = 0; i < nndata;i++) {
+//fprintf(stderr , "\n pre process v[i] data extri_0 %d extra_1 %d\n",extra_0, extra_1); //for (i = 0; i < nndata;i++) {
 //fprintf(stderr," %d:%d ", i, v[i]);
 //}
 //fprintf(stderr,"\n");
@@ -732,6 +731,25 @@ int complex_grib_out(unsigned char **sec, float *data, unsigned int ndata,
 printf("2: vmx %d vmn %d nbits %d\n", vmx, vmn, find_nbits(vmx-vmn+has_undef));
 #endif
 
+    /*
+       c2: values are deltas, which are both plus and minus
+           so range of values is up to twice range of input (one more bit)
+       c3: values are deltas of the deltas
+           the range of deltas are up to 2x initial values
+           the potential rannge of deltas of deltas are 4x initial values
+     */
+
+    j = find_nbits(vmx-vmn+has_undef);
+    if (j > 25) {
+       free(v);
+       free(sec5);
+       free(sec6);
+       max_bits -= 1;
+       if (packing_mode == 3) max_bits -= 1;
+       return complex_grib_out(sec, data, ndata, use_scale, dec_scale, bin_scale, 
+          wanted_bits, max_bits, packing_mode, use_bitmap, out);
+    }
+
     if (has_undef) {
 #pragma omp parallel for private(i)
 	for (i = 0; i < nndata; i++) {