Merge pull request #1337 from lesserwhirls/gini

Correctly handle GINI negative calibration values

Author: lesserwhirls

cdm-test/src/test/java/ucar/nc2/iosp/gini/TestGiniNegativeCalibrationValues.java

@@ -0,0 +1,37 @@
+package ucar.nc2.iosp.gini;
+
+import java.io.IOException;
+import org.junit.Assert;
+import org.junit.Test;
+import org.junit.experimental.categories.Category;
+import ucar.ma2.Array;
+import ucar.ma2.MAMath;
+import ucar.ma2.MAMath.MinMax;
+import ucar.nc2.NetcdfFile;
+import ucar.nc2.Variable;
+import ucar.unidata.util.test.TestDir;
+import ucar.unidata.util.test.category.NeedsCdmUnitTest;
+
+@Category(NeedsCdmUnitTest.class)
+public class TestGiniNegativeCalibrationValues {
+
+  private static final double comparisonTolerance = 1e-6;
+
+  @Test
+  public void testMinMaxValues() throws IOException {
+    MinMax expectedMinMax = new MinMax(-30.0, 60.0);
+    try (NetcdfFile ncfile =
+        NetcdfFile.open(TestDir.cdmUnitTestDir + "formats/gini/images_sat_NEXRCOMP_1km_n0r_n0r_20200907_0740")) {
+      Variable variable = ncfile.findVariable("Reflectivity");
+      Array array = variable.read();
+      MinMax minMax = MAMath.getMinMax(array);
+      // If the bug reported in https://github.com/Unidata/netcdf-java/issues/480 shows up, and we are
+      // incorrectly reading the calibration coefficients from the GINI file headers "Unidata cal block"
+      // (in the case of the original bug, we were decoding negative numbers as intended), the data values
+      // will be on the order of -3.9 x 10^5 for the minimum, and -2.1 x 10^5 for the maximum. Those
+      // values (radar reflectivity) should be -30 dBZ to 60 dBZ.
+      Assert.assertTrue(Math.abs(minMax.min - expectedMinMax.min) < comparisonTolerance);
+      Assert.assertTrue(Math.abs(minMax.max - expectedMinMax.max) < comparisonTolerance);
+    }
+  }
+}

clcommon/src/main/java/ucar/nc2/iosp/gini/Giniheader.java

@@ -28,13 +28,19 @@
  */
 
 class Giniheader {
+  // Unidata Gini table:
+  // https://github.com/Unidata/gempak/blob/master/gempak/tables/unidata/nex2gini.tbl
+  // Special instructions for handling negative numbers
+  // https://github.com/Unidata/gempak/blob/master/gempak/source/gemlib/mv/Linux/mvitob.c
+  //
+  private static org.slf4j.Logger log = org.slf4j.LoggerFactory.getLogger(Giniheader.class);
+
   static private final int GINI_PIB_LEN = 21;   // gini product identification block
   static private final int GINI_PDB_LEN = 512;  // gini product description block
   static private final int GINI_HED_LEN = GINI_PDB_LEN + GINI_PIB_LEN;  // gini product header
   static private final double DEG_TO_RAD = 0.017453292;
   private boolean debug = false;
   private ucar.nc2.NetcdfFile ncfile;
-  static private org.slf4j.Logger log = org.slf4j.LoggerFactory.getLogger(Giniheader.class);
   int dataStart = 0; // where the data starts
   protected int Z_type = 0;
 
@@ -379,7 +385,7 @@ void read(ucar.unidata.io.RandomAccessFile raf, ucar.nc2.NetcdfFile ncfile) thro
         break;
 
       default:
-        System.out.println("unimplemented projection");
+        log.warn("GINI: Unimplemented projection!");
     }
 
     this.ncfile.addAttribute(null, new Attribute("title", gini_GetEntityID(ent_id)));
@@ -492,6 +498,7 @@ void read(ucar.unidata.io.RandomAccessFile raf, ucar.nc2.NetcdfFile ncfile) thro
     xaxis.setDataType(DataType.DOUBLE);
     xaxis.setDimensions("x");
     xaxis.addAttribute(new Attribute(CDM.LONG_NAME, "projection x coordinate"));
+    xaxis.addAttribute(new Attribute(CF.STANDARD_NAME, "projection_x_coordinate"));
     xaxis.addAttribute(new Attribute(CDM.UNITS, "km"));
     xaxis.addAttribute(new Attribute(_Coordinate.AxisType, "GeoX"));
     double[] data = new double[nx];
@@ -520,6 +527,7 @@ void read(ucar.unidata.io.RandomAccessFile raf, ucar.nc2.NetcdfFile ncfile) thro
     yaxis.setDataType(DataType.DOUBLE);
     yaxis.setDimensions("y");
     yaxis.addAttribute(new Attribute(CDM.LONG_NAME, "projection y coordinate"));
+    yaxis.addAttribute(new Attribute(CF.STANDARD_NAME, "projection_y_coordinate"));
     yaxis.addAttribute(new Attribute(CDM.UNITS, "km"));
     yaxis.addAttribute(new Attribute(_Coordinate.AxisType, "GeoY"));
     data = new double[ny];
@@ -556,6 +564,8 @@ void read(ucar.unidata.io.RandomAccessFile raf, ucar.nc2.NetcdfFile ncfile) thro
     ncfile.addVariable(null, ct);
     ncfile.addAttribute(null, new Attribute(CDM.CONVENTIONS, _Coordinate.Convention));
 
+    var.addAttribute(new Attribute(CF.GRID_MAPPING, ct.getFullName()));
+
     // finish
     ncfile.finish();
   }
@@ -581,10 +591,35 @@ int[] getCalibrationInfo(ByteBuffer bos, int phys_elem, int ent_id) {
         for (int i = 0; i < calcod; i++) {
 
           bos.position(56 + i * 16);
+          // We expect these to always be positive, so just read the int as-is
           int minb = bos.getInt() / 10000;        /* min brightness values         */
           int maxb = bos.getInt() / 10000;       /* max brightness values         */
+
+          // Since these are data values, we need to consider the case of them being positive or negative
+          // Careful though, it's not that easy. Step one is to read the 32-bits like before...
           int mind = bos.getInt();               /* min data values               */
+          // Now for the fun.
+          // According to https://github.com/Unidata/gempak/blob/master/gempak/source/gemlib/mv/Linux/mvitob.c,
+          // which is used when writing gini files, if the first bit is set, that indicates that we have a negative
+          // value and need to do a bit more work.
+          // Check if left-most bit is 1. Right shift bit pattern of value 31 times. This will make all bits 0, or
+          // all bits 1. If all 0, then the result will equal zero. If all 1's, then the result will equal -1, and
+          // that's how we know we need to do more work.
+          if ((mind >> 31) == -1) {
+            // Set the first bit to zero, and negate the value (again, described in mvitob.c from gempak)
+            // 0x7FFFFFFF -> left-most bit is zero, all the rest are 1's.
+            // The bit-wise & results in flipping the first bit of "mind" (because we know it is 1 at this point, and
+            // 1 & 0 -> 0) while retaining the rest of the pattern of mind (because 0 & 1 -> 0, 1 & 1 -> 1).
+            // To negate, just use the negative sign...we could do ~(mind & 0x7FFFFFFF) + 1, but let's not hang out
+            // in bit operator land longer than we need to.
+            mind = -(mind & 0x7FFFFFFF);
+          }
+
+          // same thing as above
           int maxd = bos.getInt();               /* max data values               */
+          if ((maxd >> 31) == -1) {
+            maxd = -(maxd & 0x7FFFFFFF);
+          }
 
           int idscal = 1;
           while (mind % idscal == 0 && maxd % idscal == 0) {