diff --git a/mapies/mapies.py b/mapies/mapies.py
index 7bde689c25f5c3df00f1de3fe2e7556fe896759f..a944162c9f945e781b23be0293003038a1934cea 100644
--- a/mapies/mapies.py
+++ b/mapies/mapies.py
@@ -238,13 +238,13 @@ class MAPIES:
 
 
             if save:
-                if save_format == "flatten":
+                if save_format == "flattened":
                     mid_time_array = generate_mid_time_array(self.start_date, self.end_date, self.obs)
 
-                    self.to_netCDF(self.obs, self.lat, self.lon, mid_time_array, self.start_date, self.end_date)
+                    self.to_netCDF(self.obs, self.lat, self.lon, self.count_obs, mid_time_array, self.start_date, self.end_date)
                     logger.debug(f'Saved processed data to netCDF')
                 else:
-                    logger.error(f"Error: save_format '{save_format}' is not supported. Use 'flatten' instead.")
+                    logger.error(f"Error: save_format '{save_format}' is not supported. Use 'flattened' instead.")
                     raise ValueError(f"Unsupported save_format: {save_format}")
                 
 
@@ -456,6 +456,7 @@ class MAPIES:
         obs: NDArray | None = None, 
         lat: NDArray | None = None, 
         lon: NDArray | None = None, 
+        num_obs: NDArray | None = None,
         time: NDArray | None = None, 
         start_time: pd.Timestamp | None = None, 
         end_time: pd.Timestamp | None = None,
@@ -482,7 +483,9 @@ class MAPIES:
                 "lon": (["time"], lon, VARS_ATTRIBUTES['lon'].copy()),
                 "lat": (["time"], lat, VARS_ATTRIBUTES['lat'].copy()),
                 VARS_ATTRIBUTES[self.obs_var]["mapies_variable"]: (["time"], obs, VARS_ATTRIBUTES[self.obs_var].copy()),
+                "num_obs": (["time"], num_obs, VARS_ATTRIBUTES['num_obs'].copy()),
             },
+
             attrs={
                 "title": f"{VARS_ATTRIBUTES['title']['description']} from {start_time} to {end_time}",
                 "institution": "Barcelona Supercomputing Center",
@@ -500,6 +503,7 @@ class MAPIES:
             "lon": {"dtype": "float32"},
             "lat": {"dtype": "float32"},
             VARS_ATTRIBUTES[self.obs_var]["mapies_variable"]: {"dtype": "float32"}, 
+            "num_obs": {"dtype": "int32"},
         }
 
         # Save to NetCDF file
diff --git a/mapies/util/func_tools.py b/mapies/util/func_tools.py
index ff84170c4113b43b828b5e18255543024319f776..1fb6fc58e753889ce37c1687ced6ac7aba0e64e0 100644
--- a/mapies/util/func_tools.py
+++ b/mapies/util/func_tools.py
@@ -586,13 +586,17 @@ def to_netCDF_gridded_regular(
     """
 
     obs_grid = obs.astype("float32")
+    count_grid = count.astype("int32")
     obs_grid[count == 0] = np.nan
     obs_grid = np.ma.masked_invalid(obs_grid)
 
     fill_value = -9999.0
     obs_filled = np.ma.filled(obs_grid, fill_value)
+    count_grid[count==0] = fill_value
+    count_filled = np.ma.filled(count_grid, fill_value)
 
     obs_filled = obs_filled[np.newaxis,:,:]
+    count_filled = count_filled[np.newaxis,:,:]
 
     grid_lat2d = lat.astype("float32")
     grid_lon2d = lon.astype("float32")
@@ -625,6 +629,15 @@ def to_netCDF_gridded_regular(
             "_FillValue": fill_value,
             "missing_value": fill_value,
         }),
+            "num_obs": (["time", "lat", "lon"], count_filled, {
+            "units": "1", 
+            "long_name": "Number of observations used to compute the AOD550",
+            "standard_name": "Number of observations used to compute the AOD550",
+            "coordinates": "Time Latitude Longitude",
+            "_FillValue": fill_value,
+            "missing_value": fill_value,
+        }),
+
         },
         attrs={
             "title": f"aod550 averaged around {date}",
@@ -664,9 +677,14 @@ def to_netCDF_gridded_rotated(
     fill_values = -9999.0
     obs_filled = np.ma.filled(obs_grid, fill_values)
 
+    count_grid = count.astype("int32")
+    count_grid[count==0] = fill_values
+    count_filled = np.ma.filled(count_grid, fill_values)
+
     t = pd.to_datetime(date, format="%Y%m%d%H%M")
 
     obs_filled = obs_filled[np.newaxis,:,:]
+    count_filled = count_filled[np.newaxis,:,:]
 
     ds = xr.Dataset(
         coords={
@@ -708,7 +726,16 @@ def to_netCDF_gridded_rotated(
                 "grid_mapping_name": "rotated_latitude_longitude",
                 "grid_north_pole_latitude": 90-centre_lat,
                 "grid_north_pole_longitude": -180+centre_lon,
-            })
+            }),
+            "num_obs": (["time", "rlat", "rlon"], count_filled, {
+                "grid_mapping": "rotated_pole",
+                "units": "1", 
+                "long_name": "Number of observations used to compute the AOD550",
+                "standard_name": "Number of observations used to compute the AOD550",
+                "coordinates": "Time lat lon",
+                "_FillValue": fill_values,
+                "missing_value": fill_values,
+            }),
         },
         attrs={"title": f"aod550 averaged around {date}",
             "institution": "Barcelona Supercomputing Center",
diff --git a/mapies/util/variables_names.py b/mapies/util/variables_names.py
index 7e79878eb40bacc8c2d46a7b6d1216d05d5883b6..fb30d5bc77e700516798719f27b0d689ab0b8fcf 100644
--- a/mapies/util/variables_names.py
+++ b/mapies/util/variables_names.py
@@ -2,6 +2,7 @@ VARS_ATTRIBUTES={
     "lon": {"description": "Longitude of observations", "units": "degrees_east"},
     "lat": {"description": "Latitude of observations", "units": "degrees_north"},
     "Aerosol_Optical_Thickness_550_Land_Ocean_Best_Estimate": {"description": "Aerosol_Optical_Thickness_550_Land_Ocean_Best_Estimate values", "units": "unitless", "mapies_variable":"aod550"},
+    "num_obs": {"description": "Number of observations used to compute the average", "units": "unitless", "mapies_variable":"num_obs"},
     "nitrogendioxide_tropospheric_column": {"description": "nitrogendioxide_tropospheric_column", "units": "molm-2", "mapies_variable":"no2"},
     "formaldehyde_tropospheric_vertical_column": {"description": "formaldehyde_tropospheric_vertical_column", "units": "molm-2", "mapies_variable":"hcho"},
     "aerosol_index_354_388": {"description": "aerosol_index_354_388", "units": "1", "mapies_variable":"aer_ai"},
diff --git a/run/run_viirs.py b/run/run_viirs.py
index 3bfbea609051e9f4f93218fbe29b88a5e15c7b86..ed9eb22c167d71ff9e1587c89012f082cf41ca63 100644
--- a/run/run_viirs.py
+++ b/run/run_viirs.py
@@ -7,11 +7,11 @@ if __name__ == "__main__":
     # Change these variables
     # 00:30 to 02:36 from 2024-01-01 to 2024-01-03
     start_date = "202401020000"
-    end_date = "202401042359"
+    end_date = "202401022359"
 
 
 
-    outdir="/home/cgile/Documents/mapies/figures/averages/11062025"
+    outdir="/home/cgile/Documents/mapies/figures/averages/23062025"
     indir="/home/cgile/Documents/mapies/VIIRS"
 
 
@@ -28,7 +28,7 @@ if __name__ == "__main__":
     min_num_obs = 0
 
 
-    c.process_avg_data(frequency_avg=frequency_avg, min_num_obs = min_num_obs, batch_size = 100, save=True, save_format = "gridded", apply_qa = False, geospatial_crop = geospatial_crop)
+    c.process_avg_data(frequency_avg=frequency_avg, min_num_obs = min_num_obs, batch_size = 100, save=False, save_format = "gridded", apply_qa = False, geospatial_crop = geospatial_crop)
     # c.plot_2D_observations(months=[0], filename="plot_10062025.png", outdir=outdir)
     # c.plot_2D_observations(months=[1], outdir=outdir)