Correctly handle GINI negative calibration values

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) {