zaxis_compare: check bounds (bug fix)

ChangeLog

@@ -2,6 +2,10 @@
 
 	* Version 2.0.5 released
 
+2022-03-11  Uwe Schulzweida
+
+	* zaxis_compare: check bounds (bug fix)
+
 2022-03-10  Uwe Schulzweida
 
 	* ecCodes encode: fix problem with startStep for step type MIN/MAX  [report: Johannes Schick]

src/varscan.c

@@ -776,8 +776,8 @@ varDefZAxisReference(int nhlev, int nvgrid, unsigned char uuid[CDI_UUID_SIZE])
 }
 
 bool
-zaxisCompare(int zaxisID, int zaxistype, int nlevels, bool hasBounds, const double *levels, const char *longname, const char *units,
-             int ltype1, int ltype2)
+zaxis_compare(int zaxisID, int zaxistype, int nlevels, const double *levels, const double *lbounds,
+              const double *ubounds, const char *longname, const char *units, int ltype1, int ltype2)
 {
   bool differ = true;
 
@@ -786,12 +786,11 @@ zaxisCompare(int zaxisID, int zaxistype, int nlevels, bool hasBounds, const doub
   cdiInqKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_TYPEOFSECONDFIXEDSURFACE, &ltype2_0);
   const bool ltype1IsEqual = (ltype1 == ltype1_0);
   const bool ltype2IsEqual = (ltype2 == ltype2_0);
+  const bool hasBounds = (lbounds && ubounds);
 
-  printf("ltype1IsEqual, ltype2IsEqual %d %d\n", ltype1IsEqual, ltype2IsEqual);
   if (ltype1IsEqual && ltype2IsEqual && (zaxistype == zaxisInqType(zaxisID) || zaxistype == ZAXIS_GENERIC))
     {
-      const bool hasBoundsZ = (zaxisInqLbounds(zaxisID, NULL) > 0);
-      printf("   hasBoundsZ, hasBounds %d %d\n", hasBoundsZ, hasBounds);
+      const bool hasBoundsZ = (zaxisInqLbounds(zaxisID, NULL) > 0 && zaxisInqUbounds(zaxisID, NULL) > 0);
       if (nlevels == zaxisInqSize(zaxisID) && hasBoundsZ == hasBounds)
         {
           const double *dlevels = zaxisInqLevelsPtr(zaxisID);
@@ -805,6 +804,23 @@ zaxisCompare(int zaxisID, int zaxistype, int nlevels, bool hasBounds, const doub
               if (levelID == nlevels) differ = false;
             }
 
+          if (!differ && hasBounds)
+            {
+              double *bounds = malloc(2 * nlevels * sizeof(double));
+              zaxisInqLbounds(zaxisID, bounds);
+              zaxisInqUbounds(zaxisID, bounds + nlevels);
+              for (int levelID = 0; levelID < nlevels; levelID++)
+                {
+                  if (fabs(lbounds[levelID] - bounds[levelID]) > 1.e-9
+                   || fabs(ubounds[levelID] - bounds[levelID + nlevels]) > 1.e-9)
+                    {
+                      differ = true;
+                      break;
+                    }
+                }
+              free(bounds);
+            }
+
           if (!differ)
             {
               if (longname && longname[0])
@@ -833,8 +849,9 @@ struct varDefZAxisSearchState
   int resIDValue;
   int zaxistype;
   int nlevels;
-  bool hasBounds;
   const double *levels;
+  const double *lbounds;
+  const double *ubounds;
   const char *longname;
   const char *units;
   int ltype1;
@@ -846,8 +863,8 @@ varDefZAxisSearch(int id, void *res, void *data)
 {
   struct varDefZAxisSearchState *state = (struct varDefZAxisSearchState *) data;
   (void) res;
-  if (zaxisCompare(id, state->zaxistype, state->nlevels, state->hasBounds, state->levels, state->longname, state->units,
-                   state->ltype1, state->ltype2)
+  if (zaxis_compare(id, state->zaxistype, state->nlevels, state->levels, state->lbounds, state->ubounds,
+                    state->longname, state->units, state->ltype1, state->ltype2)
       == false)
     {
       state->resIDValue = id;
@@ -879,7 +896,7 @@ varDefZaxis(int vlistID, int zaxistype, int nlevels, const double *levels, const
       {
         zaxisID = vlistptr->zaxisIDs[index];
 
-        if (!zaxisCompare(zaxisID, zaxistype, nlevels, hasBounds, levels, longname, units, ltype1, ltype2))
+        if (!zaxis_compare(zaxisID, zaxistype, nlevels, levels, levels1, levels2, longname, units, ltype1, ltype2))
           {
             zaxisdefined = true;
             break;
@@ -892,7 +909,8 @@ varDefZaxis(int vlistID, int zaxistype, int nlevels, const double *levels, const
       query.zaxistype = zaxistype;
       query.nlevels = nlevels;
       query.levels = levels;
-      query.hasBounds = hasBounds;
+      query.lbounds = levels1;
+      query.ubounds = levels2;
       query.longname = longname;
       query.units = units;
       query.ltype1 = ltype1;

src/varscan.h

@@ -34,8 +34,8 @@ void varDefOptGribInt(int varID, int tile_index, long lval, const char *keyword)
 void varDefOptGribDbl(int varID, int tile_index, double dval, const char *keyword);
 int varOptGribNentries(int varID);
 
-bool zaxisCompare(int zaxisID, int zaxistype, int nlevels, bool lbounds, const double *levels, const char *longname,
-                  const char *units, int ltype1, int ltype2);
+bool zaxis_compare(int zaxisID, int zaxistype, int nlevels, const double *levels, const double *lbounds,
+                   const double *ubounds, const char *longname, const char *units, int ltype1, int ltype2);
 
 #endif
 /*

src/vlist.c

@@ -415,8 +415,9 @@ struct vgzSearchState
   int resIDValue;
   int zaxistype;
   int nlevels;
-  bool lbounds;
   const double *levels;
+  const double *lbounds;
+  const double *ubounds;
 };
 
 static enum cdiApplyRet
@@ -424,7 +425,7 @@ vgzZAxisSearch(int id, void *res, void *data)
 {
   struct vgzSearchState *state = (struct vgzSearchState *) data;
   (void) res;
-  if (zaxisCompare(id, state->zaxistype, state->nlevels, state->lbounds, state->levels, NULL, NULL, 0, -1) == false)
+  if (zaxis_compare(id, state->zaxistype, state->nlevels, state->levels, state->lbounds, state->ubounds, NULL, NULL, 0, -1) == false)
     {
       state->resIDValue = id;
       return CDI_APPLY_STOP;
@@ -449,7 +450,7 @@ vlist_generate_zaxis(int vlistID, int zaxistype, int nlevels, const double *leve
     {
       zaxisID = vlistptr->zaxisIDs[index];
 
-      if (zaxisCompare(zaxisID, zaxistype, nlevels, hasBounds, levels, NULL, NULL, 0, -1) == false)
+      if (zaxis_compare(zaxisID, zaxistype, nlevels, levels, lbounds, ubounds, NULL, NULL, 0, -1) == false)
         {
           zaxisdefined = true;
           break;
@@ -462,7 +463,6 @@ vlist_generate_zaxis(int vlistID, int zaxistype, int nlevels, const double *leve
       query.zaxistype = zaxistype;
       query.nlevels = nlevels;
       query.levels = levels;
-      query.lbounds = hasBounds;
 
       if ((zaxisglobdefined = (cdiResHFilterApply(getZaxisOps(), vgzZAxisSearch, &query) == CDI_APPLY_STOP)))
         zaxisID = query.resIDValue;