diff --git a/examples/bag_depth_has_uncertainty.py b/examples/bag_depth_has_uncertainty.py
new file mode 100644
index 0000000..38d55e3
--- /dev/null
+++ b/examples/bag_depth_has_uncertainty.py
@@ -0,0 +1,50 @@
+import os
+import logging
+
+from hyo2.abc2.lib.logging import set_logging
+from hyo2.abc2.lib.testing import Testing
+from hyo2.abc2.lib.gdal_aux import GdalAux
+from hyo2.bag.bag import BAGFile
+from hyo2.bag.meta import Meta
+from hyo2.qc4.lib.common.writers.s57_writer import S57Writer
+from hyo2.qc4.lib.common.writers.kml_writer import KmlWriter
+from hyo2.qc4.lib.common.writers.shp_writer import ShpWriter
+
+set_logging(ns_list=['hyo2.bag'])
+logger = logging.getLogger(__name__)
+
+# file_bag_0 = os.path.join(Helper.samples_folder(), "bdb_01.bag")
+file_bag_0 = r"D:\google_drive\_ccom\QC Tools\data\survey\QC Tools 4\BAG_Checks\depth_without_uncertainty\depth_no_unc.bag"
+
+root_folder = os.path.abspath(os.path.join(os.path.dirname(__file__), os.pardir))
+testing = Testing(root_folder=root_folder)
+
+GdalAux.push_gdal_error_handler()
+GdalAux.check_gdal_data(verbose=True)
+
+if os.path.exists(file_bag_0):
+    logger.debug("- file_bag_0: %s" % file_bag_0)
+
+bag_0 = BAGFile(file_bag_0, mode="r")
+logger.debug(bag_0)
+
+meta = Meta(bag_0.metadata())
+
+ret = bag_0.depth_has_uncertainty()
+logger.debug(ret)
+
+if not ret:
+    logger.info("No flags")
+    exit()
+
+s57_bn_path = os.path.join(testing.output_data_folder(), "%s.blue_notes.000" % os.path.basename(file_bag_0))
+s57_ss_path = os.path.join(testing.output_data_folder(), "%s.soundings.000" % os.path.basename(file_bag_0))
+
+flags_for_blue_notes = list()
+for entry in ret:
+    flags_for_blue_notes.append([entry[0], entry[1], "%.2f" % entry[2]])
+
+S57Writer.write_bluenotes(feature_list=flags_for_blue_notes, path=s57_bn_path, list_of_list=False)
+S57Writer.write_soundings(feature_list=ret, path=s57_ss_path, list_of_list=False)
+KmlWriter().write_bluenotes(feature_list=ret, path=s57_ss_path[:-4], list_of_list=False)
+ShpWriter().write_bluenotes(feature_list=ret, path=s57_ss_path[:-4], list_of_list=False)
diff --git a/examples/bag_uncertainty_has_depth.py b/examples/bag_uncertainty_has_depth.py
new file mode 100644
index 0000000..8673e14
--- /dev/null
+++ b/examples/bag_uncertainty_has_depth.py
@@ -0,0 +1,50 @@
+import os
+import logging
+
+from hyo2.abc2.lib.logging import set_logging
+from hyo2.abc2.lib.testing import Testing
+from hyo2.abc2.lib.gdal_aux import GdalAux
+from hyo2.bag.bag import BAGFile
+from hyo2.bag.meta import Meta
+from hyo2.qc4.lib.common.writers.s57_writer import S57Writer
+from hyo2.qc4.lib.common.writers.kml_writer import KmlWriter
+from hyo2.qc4.lib.common.writers.shp_writer import ShpWriter
+
+set_logging(ns_list=['hyo2.bag'])
+logger = logging.getLogger(__name__)
+
+# file_bag_0 = os.path.join(Helper.samples_folder(), "bdb_01.bag")
+file_bag_0 = r"D:\google_drive\_ccom\QC Tools\data\survey\QC Tools 4\BAG_Checks\uncertainty_without_depth\unc_no_depth.bag"
+
+root_folder = os.path.abspath(os.path.join(os.path.dirname(__file__), os.pardir))
+testing = Testing(root_folder=root_folder)
+
+GdalAux.push_gdal_error_handler()
+GdalAux.check_gdal_data(verbose=True)
+
+if os.path.exists(file_bag_0):
+    logger.debug("- file_bag_0: %s" % file_bag_0)
+
+bag_0 = BAGFile(file_bag_0, mode="r")
+logger.debug(bag_0)
+
+meta = Meta(bag_0.metadata())
+
+ret = bag_0.uncertainty_has_depth()
+logger.debug(ret)
+
+if not ret:
+    logger.info("No flags")
+    exit()
+
+s57_bn_path = os.path.join(testing.output_data_folder(), "%s.blue_notes.000" % os.path.basename(file_bag_0))
+s57_ss_path = os.path.join(testing.output_data_folder(), "%s.soundings.000" % os.path.basename(file_bag_0))
+
+flags_for_blue_notes = list()
+for entry in ret:
+    flags_for_blue_notes.append([entry[0], entry[1], "%.2f" % entry[2]])
+
+S57Writer.write_bluenotes(feature_list=flags_for_blue_notes, path=s57_bn_path, list_of_list=False)
+S57Writer.write_soundings(feature_list=ret, path=s57_ss_path, list_of_list=False)
+KmlWriter().write_bluenotes(feature_list=ret, path=s57_ss_path[:-4], list_of_list=False)
+ShpWriter().write_bluenotes(feature_list=ret, path=s57_ss_path[:-4], list_of_list=False)
diff --git a/examples/bag_vr_depth_has_uncertainty.py b/examples/bag_vr_depth_has_uncertainty.py
new file mode 100644
index 0000000..97596fc
--- /dev/null
+++ b/examples/bag_vr_depth_has_uncertainty.py
@@ -0,0 +1,50 @@
+import os
+import logging
+
+from hyo2.abc2.lib.logging import set_logging
+from hyo2.abc2.lib.testing import Testing
+from hyo2.abc2.lib.gdal_aux import GdalAux
+from hyo2.bag.bag import BAGFile
+from hyo2.bag.meta import Meta
+from hyo2.qc4.lib.common.writers.s57_writer import S57Writer
+from hyo2.qc4.lib.common.writers.kml_writer import KmlWriter
+from hyo2.qc4.lib.common.writers.shp_writer import ShpWriter
+
+set_logging(ns_list=['hyo2.bag'])
+logger = logging.getLogger(__name__)
+
+# file_bag_0 = os.path.join(Helper.samples_folder(), "bdb_01.bag")
+file_bag_0 = r"D:\google_drive\_ccom\QC Tools\data\survey\BAG Checks\H13275_MB_VR_MLLW.bag"
+
+root_folder = os.path.abspath(os.path.join(os.path.dirname(__file__), os.pardir))
+testing = Testing(root_folder=root_folder)
+
+GdalAux.push_gdal_error_handler()
+GdalAux.check_gdal_data(verbose=True)
+
+if os.path.exists(file_bag_0):
+    logger.debug("- file_bag_0: %s" % file_bag_0)
+
+bag_0 = BAGFile(file_bag_0, mode="r")
+logger.debug(bag_0)
+
+meta = Meta(bag_0.metadata())
+
+ret = bag_0.vr_depth_has_uncertainty()
+logger.debug(ret)
+
+if not ret:
+    logger.info("No flags")
+    exit()
+
+s57_bn_path = os.path.join(testing.output_data_folder(), "%s.blue_notes.000" % os.path.basename(file_bag_0))
+s57_ss_path = os.path.join(testing.output_data_folder(), "%s.soundings.000" % os.path.basename(file_bag_0))
+
+flags_for_blue_notes = list()
+for entry in ret:
+    flags_for_blue_notes.append([entry[0], entry[1], "%.2f" % entry[2]])
+
+S57Writer.write_bluenotes(feature_list=flags_for_blue_notes, path=s57_bn_path, list_of_list=False)
+S57Writer.write_soundings(feature_list=ret, path=s57_ss_path, list_of_list=False)
+KmlWriter().write_bluenotes(feature_list=ret, path=s57_ss_path[:-4], list_of_list=False)
+ShpWriter().write_bluenotes(feature_list=ret, path=s57_ss_path[:-4], list_of_list=False)
diff --git a/examples/bag_vr_uncertainty_has_depth.py b/examples/bag_vr_uncertainty_has_depth.py
new file mode 100644
index 0000000..a7ad0e9
--- /dev/null
+++ b/examples/bag_vr_uncertainty_has_depth.py
@@ -0,0 +1,50 @@
+import os
+import logging
+
+from hyo2.abc2.lib.logging import set_logging
+from hyo2.abc2.lib.testing import Testing
+from hyo2.abc2.lib.gdal_aux import GdalAux
+from hyo2.bag.bag import BAGFile
+from hyo2.bag.meta import Meta
+from hyo2.qc4.lib.common.writers.s57_writer import S57Writer
+from hyo2.qc4.lib.common.writers.kml_writer import KmlWriter
+from hyo2.qc4.lib.common.writers.shp_writer import ShpWriter
+
+set_logging(ns_list=['hyo2.bag'])
+logger = logging.getLogger(__name__)
+
+# file_bag_0 = os.path.join(Helper.samples_folder(), "bdb_01.bag")
+file_bag_0 = r"D:\google_drive\_ccom\QC Tools\data\survey\BAG Checks\H13275_MB_VR_MLLW.bag"
+
+root_folder = os.path.abspath(os.path.join(os.path.dirname(__file__), os.pardir))
+testing = Testing(root_folder=root_folder)
+
+GdalAux.push_gdal_error_handler()
+GdalAux.check_gdal_data(verbose=True)
+
+if os.path.exists(file_bag_0):
+    logger.debug("- file_bag_0: %s" % file_bag_0)
+
+bag_0 = BAGFile(file_bag_0, mode="r")
+logger.debug(bag_0)
+
+meta = Meta(bag_0.metadata())
+
+ret = bag_0.vr_uncertainty_has_depth()
+logger.debug(ret)
+
+if not ret:
+    logger.info("No flags")
+    exit()
+
+s57_bn_path = os.path.join(testing.output_data_folder(), "%s.blue_notes.000" % os.path.basename(file_bag_0))
+s57_ss_path = os.path.join(testing.output_data_folder(), "%s.soundings.000" % os.path.basename(file_bag_0))
+
+flags_for_blue_notes = list()
+for entry in ret:
+    flags_for_blue_notes.append([entry[0], entry[1], "%.2f" % entry[2]])
+
+S57Writer.write_bluenotes(feature_list=flags_for_blue_notes, path=s57_bn_path, list_of_list=False)
+S57Writer.write_soundings(feature_list=ret, path=s57_ss_path, list_of_list=False)
+KmlWriter().write_bluenotes(feature_list=ret, path=s57_ss_path[:-4], list_of_list=False)
+ShpWriter().write_bluenotes(feature_list=ret, path=s57_ss_path[:-4], list_of_list=False)
diff --git a/hyo2/bag/__init__.py b/hyo2/bag/__init__.py
index c24c9e0..eba1062 100644
--- a/hyo2/bag/__init__.py
+++ b/hyo2/bag/__init__.py
@@ -4,7 +4,7 @@
 """
 
 name = 'BAG'
-__version__ = '1.2.6'
+__version__ = '1.2.7'
 __author__ = 'gmasetti@ccom.unh.edu'
 __license__ = 'LGPLv3 license'
 __copyright__ = 'Copyright (c) 2024, University of New Hampshire, Center for Coastal and Ocean Mapping'
diff --git a/hyo2/bag/bag.py b/hyo2/bag/bag.py
index 2d50d0f..024638b 100644
--- a/hyo2/bag/bag.py
+++ b/hyo2/bag/bag.py
@@ -328,6 +328,110 @@ def uncertainty_greater_than(self, th: float) -> list[list[int | float]]:
 
         return xyz
 
+    def uncertainty_has_depth(self) -> list[list[int | float]]:
+        rows, cols = self.uncertainty_shape()
+        # logger.debug('shape: %s, %s' % (rows, cols))
+
+        self.populate_metadata()
+
+        x_min = self.meta.sw[0]
+        y_min = self.meta.sw[1]
+        x_res = self.meta.res_x
+        y_res = self.meta.res_y
+        # logger.debug("info: %f %f %f %f" % (x_min, y_min, x_res, y_res))
+
+        in_srs = osr.SpatialReference()
+        in_srs.ImportFromWkt(self.meta.wkt_srs)
+        out_srs = osr.SpatialReference()
+        out_srs.ImportFromEPSG(4326)
+        out_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+        ctr = osr.CoordinateTransformation(in_srs, out_srs)
+
+        mem_row = cols * 32 / 1024 / 1024
+        # mem = mem_row * rows
+        # logger.debug('estimated memory: %.1f MB' % mem)
+        chunk_size = 8096
+        chunk_rows = int(chunk_size / mem_row) + 1
+        # logger.debug('nr of rows per chunk: %s' % chunk_rows)
+
+        xyz = list()
+        for start in range(0, rows, chunk_rows):
+            stop = start + chunk_rows
+            if stop > rows:
+                stop = rows
+
+            unc = self.uncertainty(row_range=slice(start, stop))
+            # logger.info(unc)
+            dep = self.elevation(row_range=slice(start, stop))
+            # logger.info(dep)
+            unc[np.isfinite(dep)] = np.nan
+            # logger.info(unc)
+            ijs = np.argwhere(np.isfinite(unc))
+            # logger.info(ijs)
+            for ij in ijs:
+                i = ij[0]
+                j = ij[1]
+                e = x_min + j * x_res
+                n = y_min + (start + i) * y_res
+                lat, lon, _ = ctr.TransformPoint(e, n)
+                u = float(unc[i, j])
+                xyz.append([float(lat), float(lon), u])
+                # logger.info("%d,%d: %.7f %.7f %.3f" % ((start + i), j, xyz[-1][0], xyz[-1][1], xyz[-1][2]))
+
+        return xyz
+
+    def depth_has_uncertainty(self) -> list[list[int | float]]:
+        rows, cols = self.uncertainty_shape()
+        # logger.debug('shape: %s, %s' % (rows, cols))
+
+        self.populate_metadata()
+
+        x_min = self.meta.sw[0]
+        y_min = self.meta.sw[1]
+        x_res = self.meta.res_x
+        y_res = self.meta.res_y
+        # logger.debug("info: %f %f %f %f" % (x_min, y_min, x_res, y_res))
+
+        in_srs = osr.SpatialReference()
+        in_srs.ImportFromWkt(self.meta.wkt_srs)
+        out_srs = osr.SpatialReference()
+        out_srs.ImportFromEPSG(4326)
+        out_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+        ctr = osr.CoordinateTransformation(in_srs, out_srs)
+
+        mem_row = cols * 32 / 1024 / 1024
+        # mem = mem_row * rows
+        # logger.debug('estimated memory: %.1f MB' % mem)
+        chunk_size = 8096
+        chunk_rows = int(chunk_size / mem_row) + 1
+        # logger.debug('nr of rows per chunk: %s' % chunk_rows)
+
+        xyz = list()
+        for start in range(0, rows, chunk_rows):
+            stop = start + chunk_rows
+            if stop > rows:
+                stop = rows
+
+            unc = self.uncertainty(row_range=slice(start, stop))
+            # logger.info(unc)
+            dep = self.elevation(row_range=slice(start, stop))
+            # logger.info(dep)
+            dep[np.isfinite(unc)] = np.nan
+            # logger.info(dep)
+            ijs = np.argwhere(np.isfinite(dep))
+            # logger.info(ijs)
+            for ij in ijs:
+                i = ij[0]
+                j = ij[1]
+                e = x_min + j * x_res
+                n = y_min + (start + i) * y_res
+                lat, lon, _ = ctr.TransformPoint(e, n)
+                z = -float(dep[i, j])
+                xyz.append([float(lat), float(lon), z])
+                # logger.info("%d,%d: %.7f %.7f %.3f" % ((start + i), j, xyz[-1][0], xyz[-1][1], xyz[-1][2]))
+
+        return xyz
+
     def vr_uncertainty_min_max(self) -> Tuple[float, float]:
         # rows, cols = self.vr_refinements_shape()
         # logger.debug('shape: %s, %s' % (rows, cols))
@@ -391,7 +495,129 @@ def vr_uncertainty_greater_than(self, th: float) -> list[list[int | float]]:
                     # logger.debug("%d > %d,%d" % (ir_idx, rfn_r, rfn_c))
                     e = x_min + (sg_c - 0.5) * x_res + vr_ixs[sg_r, sg_c][5] + rfn_c * vr_ixs[sg_r, sg_c][3]
                     n = y_min + (sg_r - 0.5) * y_res + vr_ixs[sg_r, sg_c][6] + rfn_r * vr_ixs[sg_r, sg_c][4]
-                    lat,lon, _ = ctr.TransformPoint(e, n)
+                    lat, lon, _ = ctr.TransformPoint(e, n)
+                    xyz.append([float(lat), float(lon), unc])
+                i += ir
+
+        return xyz
+
+    def vr_depth_has_uncertainty(self) -> list[list[int | float]]:
+        # rows, cols = self.vr_refinements_shape()
+        # logger.debug('shape: %s, %s' % (rows, cols))
+
+        self.populate_metadata()
+
+        x_min = self.meta.sw[0]
+        y_min = self.meta.sw[1]
+        x_res = self.meta.res_x
+        y_res = self.meta.res_y
+        # logger.debug("info: %f %f %f %f" % (x_min, y_min, x_res, y_res))
+
+        in_srs = osr.SpatialReference()
+        in_srs.ImportFromWkt(self.meta.wkt_srs)
+        out_srs = osr.SpatialReference()
+        out_srs.ImportFromEPSG(4326)
+        out_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+        ctr = osr.CoordinateTransformation(in_srs, out_srs)
+
+        vr_unc = self[BAGFile._bag_varres_refinements][0]['depth_uncrt']
+        mask = vr_unc == BAGFile.BAG_NAN
+        vr_unc[mask] = np.nan
+        vr_dep = self[BAGFile._bag_varres_refinements][0]['depth']
+        mask = vr_dep == BAGFile.BAG_NAN
+        vr_dep[mask] = np.nan
+
+        xyz_dict = dict()
+        for idx, unc in enumerate(vr_unc):
+            dep = vr_dep[idx]
+            if np.isfinite(dep) and np.isnan(unc):
+                xyz_dict[idx] = dep
+
+        # logger.info("Located %d outliers" % len(xyz_dict))
+
+        xyz = list()
+        vr_ixs = self[BAGFile._bag_varres_metadata][:]
+        rows, cols = vr_ixs.shape
+        i = 0
+        for sg_r in range(rows):
+            for sg_c in range(cols):
+                if vr_ixs[sg_r, sg_c][1] == 0:
+                    continue
+                ir = vr_ixs[sg_r, sg_c][1] * vr_ixs[sg_r, sg_c][2]
+                for ir_idx in range(ir):
+                    j = i + ir_idx
+                    if j not in xyz_dict:
+                        continue
+                    dep = float(xyz_dict[j])
+                    # logger.debug("Located outliers: %d %f in %d,%d: %s" % (j, unc, sg_r, sg_c, vr_ixs[sg_r, sg_c]))
+                    # vr_ixs[r, c]
+                    rfn_r = ir_idx // vr_ixs[sg_r, sg_c][1]
+                    rfn_c = ir_idx % vr_ixs[sg_r, sg_c][1]
+                    # logger.debug("%d > %d,%d" % (ir_idx, rfn_r, rfn_c))
+                    e = x_min + (sg_c - 0.5) * x_res + vr_ixs[sg_r, sg_c][5] + rfn_c * vr_ixs[sg_r, sg_c][3]
+                    n = y_min + (sg_r - 0.5) * y_res + vr_ixs[sg_r, sg_c][6] + rfn_r * vr_ixs[sg_r, sg_c][4]
+                    lat, lon, _ = ctr.TransformPoint(e, n)
+                    xyz.append([float(lat), float(lon), dep])
+                i += ir
+
+        return xyz
+
+    def vr_uncertainty_has_depth(self) -> list[list[int | float]]:
+        # rows, cols = self.vr_refinements_shape()
+        # logger.debug('shape: %s, %s' % (rows, cols))
+
+        self.populate_metadata()
+
+        x_min = self.meta.sw[0]
+        y_min = self.meta.sw[1]
+        x_res = self.meta.res_x
+        y_res = self.meta.res_y
+        # logger.debug("info: %f %f %f %f" % (x_min, y_min, x_res, y_res))
+
+        in_srs = osr.SpatialReference()
+        in_srs.ImportFromWkt(self.meta.wkt_srs)
+        out_srs = osr.SpatialReference()
+        out_srs.ImportFromEPSG(4326)
+        out_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+        ctr = osr.CoordinateTransformation(in_srs, out_srs)
+
+        vr_unc = self[BAGFile._bag_varres_refinements][0]['depth_uncrt']
+        mask = vr_unc == BAGFile.BAG_NAN
+        vr_unc[mask] = np.nan
+        vr_dep = self[BAGFile._bag_varres_refinements][0]['depth']
+        mask = vr_dep == BAGFile.BAG_NAN
+        vr_dep[mask] = np.nan
+
+        xyz_dict = dict()
+        for idx, unc in enumerate(vr_unc):
+            dep = vr_dep[idx]
+            if np.isfinite(unc) and np.isnan(dep):
+                xyz_dict[idx] = unc
+
+        # logger.info("Located %d outliers" % len(xyz_dict))
+
+        xyz = list()
+        vr_ixs = self[BAGFile._bag_varres_metadata][:]
+        rows, cols = vr_ixs.shape
+        i = 0
+        for sg_r in range(rows):
+            for sg_c in range(cols):
+                if vr_ixs[sg_r, sg_c][1] == 0:
+                    continue
+                ir = vr_ixs[sg_r, sg_c][1] * vr_ixs[sg_r, sg_c][2]
+                for ir_idx in range(ir):
+                    j = i + ir_idx
+                    if j not in xyz_dict:
+                        continue
+                    unc = float(xyz_dict[j])
+                    # logger.debug("Located outliers: %d %f in %d,%d: %s" % (j, unc, sg_r, sg_c, vr_ixs[sg_r, sg_c]))
+                    # vr_ixs[r, c]
+                    rfn_r = ir_idx // vr_ixs[sg_r, sg_c][1]
+                    rfn_c = ir_idx % vr_ixs[sg_r, sg_c][1]
+                    # logger.debug("%d > %d,%d" % (ir_idx, rfn_r, rfn_c))
+                    e = x_min + (sg_c - 0.5) * x_res + vr_ixs[sg_r, sg_c][5] + rfn_c * vr_ixs[sg_r, sg_c][3]
+                    n = y_min + (sg_r - 0.5) * y_res + vr_ixs[sg_r, sg_c][6] + rfn_r * vr_ixs[sg_r, sg_c][4]
+                    lat, lon, _ = ctr.TransformPoint(e, n)
                     xyz.append([float(lat), float(lon), unc])
                 i += ir