diff --git a/src/compo/viirs_aod2ioda.py b/src/compo/viirs_aod2ioda.py
index ecf9ea8bc..2330edd11 100755
--- a/src/compo/viirs_aod2ioda.py
+++ b/src/compo/viirs_aod2ioda.py
@@ -163,6 +163,8 @@ def get_noaa_data(self):
         self.errs = self.errs[valid_pts]
         self.qcfs = self.qcfs[valid_pts]
         self.lsfs = np.zeros_like(self.lats, dtype=np.int32)
+        if np.count_nonzero(valid_pts) == 0:
+            return np.count_nonzero(valid_pts)
 
         # QCPath is the flag for retrieval path. The valid range is 0-127 in the
         # ATBD: https://www.star.nesdis.noaa.gov/jpss/documents/ATBD/ATBD_EPS_Aerosol_AOD_v3.4.pdf.
@@ -180,6 +182,7 @@ def get_noaa_data(self):
             self.errs = 0.111431 + 0.128699 * self.vals    # over land (dark)
             self.errs[qcpath % 2 == 1] = 0.00784394 + 0.219923 * self.vals[qcpath % 2 == 1]  # over ocean
             self.errs[qcpath % 4 == 2] = 0.0550472 + 0.299558 * self.vals[qcpath % 4 == 2]   # over bright land
+        return np.count_nonzero(valid_pts)
 
     def get_nasa_dt_data(self):
         # For NASA Dark Target
@@ -206,6 +209,7 @@ def get_nasa_dt_data(self):
 
         # Flip QC flags for PreQC (0->3, 3->0)
         self.qcfs = nasa_flip_qc[self.qcfs[valid_pts].astype(np.int32)]
+        return np.count_nonzero(valid_pts)
 
     def get_nasa_db_data(self):
         # For NASA Deep Blue
@@ -213,13 +217,18 @@ def get_nasa_db_data(self):
         self.lats = self.ncd.variables['Latitude'][:].ravel()
         # Only QC=3 pixels are retained with Aerosol_Optical_Thickness_550_Land_Ocean_Best_Estimate
         self.vals = self.ncd.variables['Aerosol_Optical_Thickness_550_Land_Ocean_Best_Estimate'][:].ravel()
+        eu_land = self.ncd.variables['Aerosol_Optical_Thickness_550_Expected_Uncertainty_Land'][:].ravel()
+        eu_ocean = self.ncd.variables['Aerosol_Optical_Thickness_550_Expected_Uncertainty_Ocean'][:].ravel()
 
         # Keep valid data points only
-        valid_pts = ~self.vals.mask
+        nan_mask = (np.isnan(self.vals) | np.isnan(eu_land) | np.isnan(eu_ocean))
+        valid_pts = (~self.vals.mask & (~(eu_land < 0) | ~(eu_ocean < 0))) & (~nan_mask)
         self.lons = self.lons[valid_pts]
         self.lats = self.lats[valid_pts]
         self.vals = self.vals[valid_pts]
         self.lsfs = np.zeros_like(self.lats, dtype=np.int32)
+        if np.count_nonzero(valid_pts) == 0:
+            return np.count_nonzero(valid_pts)
 
         npts_land = self.ncd.variables['Number_Of_Pixels_Used_Land'][:].ravel()
         npts_ocean = self.ncd.variables['Number_Of_Pixels_Used_Ocean'][:].ravel()
@@ -239,8 +248,9 @@ def get_nasa_db_data(self):
         # Lee et al. (2024): https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JD040082?af=R
         # PUE should fit for DA purpose better according to
         # Hsu et al. (2018): https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JD029688
-        eu_land = self.ncd.variables['Aerosol_Optical_Thickness_550_Expected_Uncertainty_Land'][:].ravel()[valid_pts]
-        eu_ocean = self.ncd.variables['Aerosol_Optical_Thickness_550_Expected_Uncertainty_Ocean'][:].ravel()[valid_pts]
+        eu_land = eu_land[valid_pts]
+        eu_ocean = eu_ocean[valid_pts]
+
         self.errs = np.ones_like(eu_land)
         qaf_land = self.ncd.variables['Aerosol_Optical_Thickness_QA_Flag_Land'][:].ravel()[valid_pts]
         qaf_ocean = self.ncd.variables['Aerosol_Optical_Thickness_QA_Flag_Ocean'][:].ravel()[valid_pts]
@@ -270,6 +280,7 @@ def get_nasa_db_data(self):
             # VIIRS DeepBlue Expected Error (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JD029688)
             self.errs = np.add(0.05, np.multiply(0.2, self.vals))
             AttrData['errorMethod'] = 'Expected Error (EE)'
+        return np.count_nonzero(valid_pts)
 
     def read(self):
         # Make empty lists for the output vars
@@ -309,8 +320,11 @@ def read(self):
             # Get the platform and sensor name
             self.get_platform_sensor_names()
 
-            # Get VIIRS data
-            get_viirs_data()
+            # Get VIIRS data and skip the file if no valid data
+            nvalid_pts = get_viirs_data()
+            if nvalid_pts == 0:
+                print(f"  No valid pixels, skip {f}")
+                continue
 
             # assign the observation time based on time coverage
             self.obs_time = np.full(np.shape(self.lons), round(0.5*(self.s_time + self.e_time)), dtype=np.int64)
@@ -401,13 +415,14 @@ def main():
 
     args = parser.parse_args()
 
-    args_in_dict = {'input': args.input,
-                    'error_method': args.error_method,
-                    'provider': args.provider,
-                    'retrieval_method': args.retrieval_method,
-                    'thin': args.thin,
-                    'date_range': args.date_range,
-                    }
+    args_in_dict = {
+        'input': args.input,
+        'error_method': args.error_method,
+        'provider': args.provider,
+        'retrieval_method': args.retrieval_method,
+        'thin': args.thin,
+        'date_range': args.date_range,
+    }
 
     # setup the IODA writer
 
diff --git a/test/testoutput/viirs_db_aod.nc b/test/testoutput/viirs_db_aod.nc
index 4a4b34ccc..a4792174b 100644
--- a/test/testoutput/viirs_db_aod.nc
+++ b/test/testoutput/viirs_db_aod.nc
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:ba5b7921736215ab589247c2682ffb66cdb25909a504ccdbed42135646c7de05
+oid sha256:26b534141f2d5ac8516f07eef94da578421b71ffb8b0343c155112dd903a8e44
 size 15256
diff --git a/test/testoutput/viirs_dt_aod.nc b/test/testoutput/viirs_dt_aod.nc
index 0573a0540..1f13187bf 100644
--- a/test/testoutput/viirs_dt_aod.nc
+++ b/test/testoutput/viirs_dt_aod.nc
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:9f20907b39f3eadba760aab94c934c9cf7ec6dd153a45ea0971be0b9cabe9080
+oid sha256:b16bd4a1882214fd2f51ebb30c3c41cf342b55b869f0ae70ade695a8ebfd195f
 size 15256