isosurface: added docu.

ChangeLog

+2020-04-09 Uwe Schulzweida
+
+	* New operator isosurface: Extract isosurface
+
 2020-03-26 Uwe Schulzweida
 
 	* New operator addtrend: Add trend

NEWS

@@ -6,9 +6,13 @@ Improvement
 Version 1.9.9 (29 October 2020):
 
    New features:
+     * showattribute: added wildcard support
+     * Diff: added option maxcount=<num>: Stop after num different fields
    New operators:
      * apply: Apply an operator on each input file
      * verifygrid: Verify grid coordinates
+     * addtrend: Add trend
+     * isosurface: Extract isosurface
    Fixed bugs:
      * intlevel3d: wrong result since v1.9.4 [Bug #9468]
      * dv2uv, uv2dv: wrong result works only on first level since v1.9.8 [Bug #9441]

OPERATORS

@@ -531,6 +531,7 @@ Operator catalog:
    Intlevel      intlevel        Linear level interpolation
    Intlevel3d    intlevel3d      Linear level interpolation onto a 3D vertical coordinate
    Intlevel3d    intlevelx3d     like intlevel3d but with extrapolation
+   Isosurface    isosurface      Extract isosurface data
    Inttime       inttime         Interpolation between timesteps
    Inttime       intntime        Interpolation between timesteps
    Intyear       intyear         Interpolation between two years

doc/tex/Modules

@@ -117,6 +117,7 @@ Vertint       Interpolation
 Vertintap     Interpolation
 Intlevel      Interpolation
 Intlevel3d    Interpolation
+Isosurface    Interpolation
 Inttime       Interpolation
 Intyear       Interpolation
 Spectral      Transformation

doc/tex/mod/Isosurface

+@BeginModule
+@NewPage
+@Name      = Isosurface
+@Title     = Extract isosurface
+@Section   = Interpolation
+@Class     = Interpolation
+@Arguments = infile outfile
+@Operators = isosurface
+
+@BeginDescription
+This operator computes isosurface data from all 3D variables at the isosurface value specified in isovalue.
+The isosurface is calculated by linear interpolation between two layers. The result is a horizonal 2D field.
+@EndDescription
+@EndModule
+
+
+@BeginOperator_isosurface
+@Title     = Extract isosurface data
+@Parameter = isovalue
+@EndOperator
+
+
+@BeginParameter
+@Item = isovalue
+FLOAT  Isosurface value
+@EndParameter

src/Isosurface.cc

@@ -66,7 +66,6 @@ void *
 Isosurface(void *process)
 {
   int nrecs;
-  int varID, levelID;
   int zaxisID1 = -1;
   size_t nmiss;
 
@@ -90,14 +89,14 @@ Isosurface(void *process)
   vlistDefTaxis(vlistID2, taxisID2);
 
   const auto nzaxis = vlistNzaxis(vlistID1);
-  int nlevel = 0;
+  int nlevels = 0;
   int i;
   for (i = 0; i < nzaxis; i++)
     {
       const auto zaxisID = vlistZaxis(vlistID1, i);
-      nlevel = zaxisInqSize(zaxisID);
+      nlevels = zaxisInqSize(zaxisID);
       if (zaxisInqType(zaxisID) != ZAXIS_HYBRID && zaxisInqType(zaxisID) != ZAXIS_HYBRID_HALF)
-        if (nlevel > 1)
+        if (nlevels > 1)
           {
             zaxisID1 = zaxisID;
             break;
@@ -105,7 +104,7 @@ Isosurface(void *process)
     }
   if (i == nzaxis) cdoAbort("No processable variable found!");
 
-  const auto nlev1 = nlevel;
+  const auto nlev1 = nlevels;
   Varray<double> lev1(nlev1);
   cdoZaxisInqLevels(zaxisID1, &lev1[0]);
 
@@ -116,7 +115,8 @@ Isosurface(void *process)
   const auto streamID2 = cdoOpenWrite(1);
   cdoDefVlist(streamID2, vlistID2);
 
-  VarList varList2;
+  VarList varList1, varList2;
+  varListInit(varList1, vlistID1);
   varListInit(varList2, vlistID2);
 
   const auto gridsizemax = vlistGridsizeMax(vlistID1);
@@ -130,11 +130,7 @@ Isosurface(void *process)
   FieldVector2D vars1;
   fieldsFromVlist(vlistID1, vars1, FIELD_VEC);
 
-  for (varID = 0; varID < nvars; varID++)
-    {
-      const auto zaxisID = vlistInqVarZaxis(vlistID1, varID);
-      liso[varID] = (zaxisID == zaxisID1);
-    }
+  for (int varID = 0; varID < nvars; varID++) liso[varID] = (varList1[varID].zaxisID == zaxisID1);
 
   int tsID = 0;
   while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
@@ -142,17 +138,18 @@ Isosurface(void *process)
       taxisCopyTimestep(taxisID2, taxisID1);
       cdoDefTimestep(streamID2, tsID);
 
-      for (varID = 0; varID < nvars; varID++) vars[varID] = false;
+      for (int varID = 0; varID < nvars; varID++) vars[varID] = false;
 
       for (int recID = 0; recID < nrecs; recID++)
         {
+          int varID, levelID;
           cdoInqRecord(streamID1, &varID, &levelID);
           cdoReadRecord(streamID1, vars1[varID][levelID].vec_d.data(), &nmiss);
           vars1[varID][0].nmiss += nmiss;
           vars[varID] = true;
         }
 
-      for (varID = 0; varID < nvars; varID++)
+      for (int varID = 0; varID < nvars; varID++)
         {
           if (vars[varID])
             {
@@ -165,7 +162,7 @@ Isosurface(void *process)
                 }
               else
                 {
-                  for (levelID = 0; levelID < varList2[varID].nlevels; levelID++)
+                  for (int levelID = 0; levelID < varList2[varID].nlevels; levelID++)
                     {
                       auto &single = vars1[varID][levelID].vec_d;
                       nmiss = varrayNumMV(varList2[varID].gridsize, single, varList2[varID].missval);

src/module_list.h

@@ -94,7 +94,7 @@ static const module_t module_Intntime       = {Intntime       , InttimeHelp
 static const module_t module_Intyear        = {Intyear        , IntyearHelp       , IntyearOperators        , EXPOSED  , CDI_REAL , 2  , OBASE , NoRestriction };
 static const module_t module_Invert         = {Invert         , InvertHelp        , InvertOperators         , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };
 static const module_t module_Invertlev      = {Invertlev      , InvertlevHelp     , InvertlevOperators      , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };
-static const module_t module_Isosurface     = {Isosurface     , {}                , IsosurfaceOperators     , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };
+static const module_t module_Isosurface     = {Isosurface     , IsosurfaceHelp    , IsosurfaceOperators     , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };
 static const module_t module_Lic            = {Lic            , {}                , LicOperators            , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };
 static const module_t module_MapReduce      = {MapReduce      , MapReduceHelp     , MapReduceOperators      , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };
 static const module_t module_Maskbox        = {Maskbox        , MaskboxHelp       , MaskboxOperators        , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };

src/operator_help.h

@@ -4304,6 +4304,21 @@ static const std::vector<std::string> Intlevel3dHelp = {
     "    infile2      STRING  target vertical coordinate field",
 };
 
+static const std::vector<std::string> IsosurfaceHelp = {
+    "NAME",
+    "    isosurface - Extract isosurface",
+    "",
+    "SYNOPSIS",
+    "    isosurface,isovalue  infile outfile",
+    "",
+    "DESCRIPTION",
+    "    This operator computes isosurface data from all 3D variables at the isosurface value specified in isovalue.",
+    "    The isosurface is calculated by linear interpolation between two layers. The result is a horizonal 2D field.",
+    "",
+    "PARAMETER",
+    "    isovalue  FLOAT  Isosurface value",
+};
+
 static const std::vector<std::string> InttimeHelp = {
     "NAME",
     "    inttime, intntime - Time interpolation",