added method to identify grid cells with high uncertainty

Author: giumas

.gitignore

@@ -62,3 +62,6 @@ docs/_build/*
 # notebooks checkpoint
 *.ipynb_checkpoints
 
+data/temp/
+data/download/
+data/output/

data/.gitkeep

@@ -0,0 +1,38 @@
+import os
+import logging
+
+from hyo2.abc2.lib.logging import set_logging
+from hyo2.abc2.lib.testing import Testing
+from hyo2.bag.bag import BAGFile
+from hyo2.bag.meta import Meta
+from hyo2.qc4.lib.common.writers.s57_writer import S57Writer
+
+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"
+th = 2.0
+
+root_folder = os.path.abspath(os.path.join(os.path.dirname(__file__), os.pardir))
+testing = Testing(root_folder=root_folder)
+
+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_greater_than(th=th)
+logger.debug(ret)
+
+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)

examples/bag_uncertainty_greater_than.py

@@ -0,0 +1,38 @@
+import os
+import logging
+
+from hyo2.abc2.lib.logging import set_logging
+from hyo2.abc2.lib.testing import Testing
+from hyo2.bag.bag import BAGFile
+from hyo2.bag.meta import Meta
+from hyo2.qc4.lib.common.writers.s57_writer import S57Writer
+
+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"
+th = 2.0
+
+root_folder = os.path.abspath(os.path.join(os.path.dirname(__file__), os.pardir))
+testing = Testing(root_folder=root_folder)
+
+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_greater_than(th=th)
+logger.debug(ret)
+
+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)

examples/bag_vr_uncertainty_greater_than.py

@@ -4,7 +4,7 @@
 """
 
 name = 'BAG'
-__version__ = '1.2.5'
+__version__ = '1.2.6'
 __author__ = '[email protected]'
 __license__ = 'LGPLv3 license'
 __copyright__ = 'Copyright (c) 2024, University of New Hampshire, Center for Coastal and Ocean Mapping'

hyo2/bag/__init__.py

@@ -2,6 +2,7 @@
 import logging
 from typing import Tuple
 
+from osgeo import osr
 import numpy as np
 from lxml import etree, isoschematron
 
@@ -86,7 +87,7 @@ def create_template(cls, name):
             elevation.attrs.create(cls._bag_elevation_min_ev, 0.0, shape=(), dtype=np.float32)
             elevation.attrs.create(cls._bag_elevation_max_ev, 0.0, shape=(), dtype=np.float32)
 
-            new_bag.create_dataset(cls._bag_metadata, shape=(1, ), dtype="S1")
+            new_bag.create_dataset(cls._bag_metadata, shape=(1,), dtype="S1")
 
             tracking_list = new_bag.create_dataset(cls._bag_tracking_list, shape=(), dtype=cls._bag_tracking_list_type)
             tracking_list.attrs.create(cls._bag_tracking_list_len, 0, shape=(), dtype=np.uint32)
@@ -281,6 +282,51 @@ def uncertainty_min_max(self) -> Tuple[float, float]:
 
         return unc_min, unc_max
 
+    def uncertainty_greater_than(self, th: float) -> 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)
+        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))
+            ijs = np.argwhere(unc > th)
+            for ij in ijs:
+                i = ij[0]
+                j = ij[1]
+                e = x_min + j * x_res
+                n = y_min + i * y_res
+                lon, lat, _ = ctr.TransformPoint(e, n)
+                u = float(unc[i, j])
+                xyz.append([float(lat), float(lon), u])
+                # logger.info("%s" % (xyz[-1]))
+
+        return xyz
+
     def vr_uncertainty_min_max(self) -> Tuple[float, float]:
         # rows, cols = self.vr_refinements_shape()
         # logger.debug('shape: %s, %s' % (rows, cols))
@@ -293,7 +339,64 @@ def vr_uncertainty_min_max(self) -> Tuple[float, float]:
 
         return np.nanmin(vr_unc), np.nanmax(vr_unc)
 
+    def vr_uncertainty_greater_than(self, th: float) -> 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)
+        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
+
+        xyz_dict = dict()
+        for idx, unc in enumerate(vr_unc):
+            if unc > th:
+                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]
+                    lon, lat, _ = ctr.TransformPoint(e, n)
+                    xyz.append([float(lat), float(lon), unc])
+                i += ir
+
+        return xyz
+
     def has_density(self):
+        # noinspection PyBroadException
         try:
             self[BAGFile._bag_elevation_solution]['num_soundings']
         except Exception:
@@ -437,7 +540,7 @@ def substitute_metadata(self, path):
 
         del self[BAGFile._bag_metadata]
         xml_sz = len(xml_string)
-        ds = self.create_dataset(self._bag_metadata, (xml_sz, ), dtype="S1")
+        ds = self.create_dataset(self._bag_metadata, (xml_sz,), dtype="S1")
         for i, bt in enumerate(xml_string):
             ds[i] = bytes([bt])
 
@@ -566,7 +669,7 @@ def modify_wkt_prj(self, wkt_hor, wkt_ver=None):
 
         new_xml = etree.tostring(xml_tree, pretty_print=True)
         del self[BAGFile._bag_metadata]
-        ds = self.create_dataset(BAGFile._bag_metadata, shape=(len(new_xml), ), dtype="S1")
+        ds = self.create_dataset(BAGFile._bag_metadata, shape=(len(new_xml),), dtype="S1")
         for i, x in enumerate(new_xml):
             ds[i] = bytes([x])