smooth: added support for missing values

src/Smooth.c

@@ -40,31 +40,14 @@ double smooth_nbr_compute_weights(unsigned num_neighbors, const int *restrict sr
 
   double dist_tot = 0.; // sum of neighbor distances (for normalizing)
 
-  if ( src_grid_mask )
+  for ( unsigned n = 0; n < num_neighbors; ++n )
     {
-      for ( unsigned n = 0; n < num_neighbors; ++n )
+      nbr_mask[n] = FALSE;
+      if ( nbr_add[n] >= 0 && src_grid_mask[nbr_add[n]] )
         {
-          nbr_mask[n] = FALSE;
-          if ( nbr_add[n] >= 0 )
-            if ( src_grid_mask[nbr_add[n]] )
-              {
-                nbr_dist[n] = 1./nbr_dist[n];
-                dist_tot += nbr_dist[n];
-                nbr_mask[n] = TRUE;
-              }
-        }
-    }
-  else
-    {
-      for ( unsigned n = 0; n < num_neighbors; ++n )
-        {
-          nbr_mask[n] = FALSE;
-          if ( nbr_add[n] >= 0 )
-            {
-              nbr_dist[n] = intlin(nbr_dist[n], weight0, 0, weightInf, search_radius);
-              dist_tot += nbr_dist[n];
-              nbr_mask[n] = TRUE;
-            }
+          nbr_dist[n] = intlin(nbr_dist[n], weight0, 0, weightInf, search_radius);
+          dist_tot += nbr_dist[n];
+          nbr_mask[n] = TRUE;
         }
     }
 
@@ -100,6 +83,13 @@ void smooth(int gridID, double missval, const double *restrict array1, double *r
   int gridID0 = gridID;
   unsigned gridsize = gridInqSize(gridID);
 
+  int *mask = (int*) Malloc(gridsize*sizeof(double));
+  for ( unsigned i = 0; i < gridsize; ++i )
+    if ( DBL_IS_EQUAL(array1[i], missval) )
+      mask[i] = 0;
+    else
+      mask[i] = 1;
+  
   double *xvals = (double*) Malloc(gridsize*sizeof(double));
   double *yvals = (double*) Malloc(gridsize*sizeof(double));
 
@@ -161,22 +151,27 @@ void smooth(int gridID, double missval, const double *restrict array1, double *r
       grid_search_nbr(gs, num_neighbors, nbr_add, nbr_dist, xvals[i], yvals[i]);
 
       /* Compute weights based on inverse distance if mask is false, eliminate those points */
-      double dist_tot = smooth_nbr_compute_weights(num_neighbors, NULL, nbr_mask, nbr_add, nbr_dist,
+      double dist_tot = smooth_nbr_compute_weights(num_neighbors, mask, nbr_mask, nbr_add, nbr_dist,
                                                    search_radius, weight0, weightInf);
 
       /* Normalize weights and store the link */
       unsigned nadds = smooth_nbr_normalize_weights(num_neighbors, dist_tot, nbr_mask, nbr_add, nbr_dist);
       if ( nadds )
         {
-          
+          /*
           printf("n %u %d nadds %u dis %g\n", i, nbr_add[0], nadds, nbr_dist[0]);
           for ( unsigned n = 0; n < nadds; ++n )
             printf("   n %u add %d dis %g\n", n, nbr_add[n], nbr_dist[n]);
-          
+          */
           double result = 0;
           for ( unsigned n = 0; n < nadds; ++n ) result += array1[nbr_add[n]]*nbr_dist[n];
           array2[i] = result;
         }
+      else
+        {
+          (*nmiss)++;
+          array2[i] = missval;
+        }
     }
 
   finish = clock();
@@ -187,6 +182,7 @@ void smooth(int gridID, double missval, const double *restrict array1, double *r
 
   if ( gridID0 != gridID ) gridDestroy(gridID);
 
+  Free(mask);
   Free(xvals);
   Free(yvals);
 }
@@ -310,8 +306,9 @@ void *Smooth(void *argument)
   int varID, levelID;
   int nmiss;
   int gridtype;
+  int nsmooth = 4;
   int max_points = 5;
-  double search_radius = 60, weight0 = 0.25, weightInf = 0.;
+  double search_radius = 120, weight0 = 0.25, weightInf = 0.25;
 
   search_radius *= DEG2RAD;
 
@@ -376,11 +373,16 @@ void *Smooth(void *argument)
 	      double missval = vlistInqVarMissval(vlistID1, varID);
 	      gridID = vlistInqVarGrid(vlistID1, varID);
 
-              if ( operatorID == SMOOTH )
-                smooth(gridID, missval, array1, array2, &nmiss, max_points, search_radius, weight0, weightInf);
-              else if ( operatorID == SMOOTH9 )
-                smooth9(gridID, missval, array1, array2, &nmiss);
-
+              for ( int i = 0; i < nsmooth; ++i )
+                {
+                  if ( operatorID == SMOOTH )
+                    smooth(gridID, missval, array1, array2, &nmiss, max_points, search_radius, weight0, weightInf);
+                  else if ( operatorID == SMOOTH9 )
+                    smooth9(gridID, missval, array1, array2, &nmiss);
+                  
+                  memcpy(array1, array2, gridsize*sizeof(double));
+                }
+          
 	      streamDefRecord(streamID2, varID, levelID);
 	      streamWriteRecord(streamID2, array2, nmiss);		
 	    }