diff --git a/utils/obsproc/Viirsaod2Ioda.h b/utils/obsproc/Viirsaod2Ioda.h
index 3bcdb26b0..086caf185 100644
--- a/utils/obsproc/Viirsaod2Ioda.h
+++ b/utils/obsproc/Viirsaod2Ioda.h
@@ -133,7 +133,8 @@ namespace gdasapp {
       std::vector<std::vector<int>> mask_s;
 
       if ( fullConfig_.has("binning") ) {
-        // deal with longitude when points cross the international date line
+        // Do superobbing
+        // Deal with longitude when points cross the international date line
         float minLon = std::numeric_limits<float>::max();
         float maxLon = std::numeric_limits<float>::min();
 
@@ -158,19 +159,17 @@ namespace gdasapp {
             }
         }
 
-
         lat2d_s = gdasapp::superobutils::subsample2D(lat, mask, fullConfig_);
         mask_s = gdasapp::superobutils::subsample2D(mask, mask, fullConfig_);
         if (fullConfig_.has("binning.cressman radius")) {
-        // Weighted-average (cressman) does not work yet. Need to develop subsample2D_cressman
-        // function to superob.h or incorporate weighted average to subsample2D function.
-        // obsvalue_s = gdasapp::superobutils::subsample2D_cressman(obsvalue, lat, lon,
-        //              lat2d_s, lon2d_s, mask, fullConfig_);
-        // obserror_s = gdasapp::superobutils::subsample2D_cressman(obserror, lat, lon,
-        // lat2d_s, lon2d_s, mask, fullConfig_);
-          obsvalue_s = gdasapp::superobutils::subsample2D(obsvalue, mask, fullConfig_);
-          obserror_s = gdasapp::superobutils::subsample2D(obserror, mask, fullConfig_);
+        // Weighted-average (cressman) superob
+          bool useCressman = true;
+          obsvalue_s = gdasapp::superobutils::subsample2D(obsvalue, mask, fullConfig_,
+                       useCressman, lat, lon, lat2d_s, lon2d_s);
+          obserror_s = gdasapp::superobutils::subsample2D(obserror, mask, fullConfig_,
+                       useCressman, lat, lon, lat2d_s, lon2d_s);
         } else {
+        // Simple-average superob
           obsvalue_s = gdasapp::superobutils::subsample2D(obsvalue, mask, fullConfig_);
           obserror_s = gdasapp::superobutils::subsample2D(obserror, mask, fullConfig_);
         }
diff --git a/utils/obsproc/superob.h b/utils/obsproc/superob.h
index 413527ea2..759014f1a 100644
--- a/utils/obsproc/superob.h
+++ b/utils/obsproc/superob.h
@@ -4,15 +4,26 @@
 #include <string>
 #include <vector>
 
+#include "atlas/util/Earth.h"
+#include "atlas/util/Geometry.h"
+#include "atlas/util/Point.h"
 #include "eckit/config/LocalConfiguration.h"
 
+
+
 namespace gdasapp {
   namespace superobutils {
     // Function to perform subsampling/binning of gridded data with a given stride
     template <typename T>
     std::vector<std::vector<T>> subsample2D(const std::vector<std::vector<T>>& inputArray,
                                             const std::vector<std::vector<int>>& mask,
-                                            const eckit::Configuration & fullConfig) {
+                                            const eckit::Configuration & fullConfig,
+                                            bool useCressman = false,
+                                            const std::vector<std::vector<T>>& inputlat = {},
+                                            const std::vector<std::vector<T>>& inputlon = {},
+                                            const std::vector<std::vector<T>>& targetlat = {},
+                                            const std::vector<std::vector<T>>& targetlon = {}
+) {
       // Get the binning configuration
       int stride;
       int minNumObs;
@@ -31,30 +42,73 @@ namespace gdasapp {
       // Perform subsampling
       T sum;
       int count;
-      for (int i = 0; i < subsampledRows; ++i) {
-        for (int j = 0; j < subsampledCols; ++j) {
-          count = 0;
-          sum = static_cast<T>(0);
-          // Compute the average within the stride
-          for (int si = 0; si < stride; ++si) {
-            for (int sj = 0; sj < stride; ++sj) {
-              int row = i * stride + si;
-              int col = j * stride + sj;
-              if (row < numRows && col < numCols && mask[row][col] == 1) {
-                sum += inputArray[row][col];
-                count++;
+      if (!useCressman) {
+        for (int i = 0; i < subsampledRows; ++i) {
+          for (int j = 0; j < subsampledCols; ++j) {
+            count = 0;
+            sum = static_cast<T>(0);
+            // Compute the average within the stride
+            for (int si = 0; si < stride; ++si) {
+              for (int sj = 0; sj < stride; ++sj) {
+                int row = i * stride + si;
+                int col = j * stride + sj;
+                if (row < numRows && col < numCols && mask[row][col] == 1) {
+                  sum += inputArray[row][col];
+                  count++;
+                }
               }
             }
+
+            // Calculate the average and store it in the subsampled array
+            if ( count < minNumObs ) {
+              subsampled[i][j] = static_cast<T>(-9999);
+            } else {
+              subsampled[i][j] = sum / static_cast<T>(count);
+            }
           }
+        }
+      } else {
+      // Apply Cressman interpolation if selected
+        // Perform Cressman interpolation
+        double cressmanRadius;
+        fullConfig.get("binning.cressman radius", cressmanRadius);
+        for (int i = 0; i < subsampledRows; ++i) {
+          for (int j = 0; j < subsampledCols; ++j) {
+            // Initialize sum and sumWeights for Cressman interpolation
+            count = 0;
+            sum = static_cast<T>(0);
+            T sumWeights = static_cast<T>(0);
+            // Loop through neighboring points for interpolation
+            for (int si = 0; si < stride; ++si) {
+              for (int sj = 0; sj < stride; ++sj) {
+                int row = i * stride + si;
+                int col = j * stride + sj;
+                if (row < numRows && col < numCols && mask[row][col] == 1) {
+                  atlas::PointLonLat p0(inputlon[row][col], inputlat[row][col]);
+                  atlas::PointLonLat p1(targetlon[i][j], targetlat[i][j]);
+                  double distance = atlas::util::Earth::distance(p0, p1)/1000.0;
 
-          // Calculate the average and store it in the subsampled array
-          if ( count < minNumObs ) {
-            subsampled[i][j] = static_cast<T>(-9999);
-          } else {
-            subsampled[i][j] = sum / static_cast<T>(count);
+                  double distance_sq = distance * distance;
+                  double cressmanRadius_sq = cressmanRadius * cressmanRadius;
+                  double weight = (distance <= cressmanRadius) ? (cressmanRadius_sq - distance_sq )
+                                  / (cressmanRadius_sq + distance_sq) : 0.0;
+                  sum += inputArray[row][col] * weight;
+                  sumWeights += weight;
+                  count++;
+                }
+              }
+            }
+
+            // Update subsampled value with Cressman interpolation
+            if (count < minNumObs || sumWeights == 0.0) {
+              subsampled[i][j] = static_cast<T>(-9999);
+            } else {
+              subsampled[i][j] = sum / static_cast<T>(sumWeights);
+            }
           }
         }
       }
+
       std::cout << " done subsampling" << std::endl;
       return subsampled;
     }