set GDAL's axes traditional order for BAG higher uncertainty

hydroffice · Jun 30, 2024 · 51ad448 · 51ad448
1 parent e47150c
commit 51ad448
Show file tree

Hide file tree

Showing 3 changed files with 25 additions and 7 deletions.
diff --git a/examples/bag_uncertainty_greater_than.py b/examples/bag_uncertainty_greater_than.py
@@ -3,9 +3,12 @@
 
 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__)
@@ -17,6 +20,9 @@
 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)
 
@@ -36,3 +42,5 @@
     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_greater_than.py b/examples/bag_vr_uncertainty_greater_than.py
@@ -3,9 +3,12 @@
 
 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__)
@@ -23,6 +26,9 @@
 bag_0 = BAGFile(file_bag_0, mode="r")
 logger.debug(bag_0)
 
+GdalAux.push_gdal_error_handler()
+GdalAux.check_gdal_data(verbose=True)
+
 meta = Meta(bag_0.metadata())
 
 ret = bag_0.vr_uncertainty_greater_than(th=th)
@@ -36,3 +42,5 @@
     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/bag.py b/hyo2/bag/bag.py
@@ -292,12 +292,13 @@ def uncertainty_greater_than(self, th: float) -> list[list[int | float]]:
         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))
+        # 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
@@ -320,7 +321,7 @@ def uncertainty_greater_than(self, th: float) -> list[list[int | float]]:
                 j = ij[1]
                 e = x_min + j * x_res
                 n = y_min + i * y_res
-                lon, lat, _ = ctr.TransformPoint(e, n)
+                lat, lon, _ = ctr.TransformPoint(e, n)
                 u = float(unc[i, j])
                 xyz.append([float(lat), float(lon), u])
                 # logger.info("%s" % (xyz[-1]))
@@ -349,12 +350,13 @@ def vr_uncertainty_greater_than(self, th: float) -> list[list[int | float]]:
         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))
+        # 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']
@@ -366,7 +368,7 @@ def vr_uncertainty_greater_than(self, th: float) -> list[list[int | float]]:
             if unc > th:
                 xyz_dict[idx] = unc
 
-        logger.info("Located %d outliers" % len(xyz_dict))
+        # logger.info("Located %d outliers" % len(xyz_dict))
 
         xyz = list()
         vr_ixs = self[BAGFile._bag_varres_metadata][:]
@@ -382,14 +384,14 @@ def vr_uncertainty_greater_than(self, th: float) -> list[list[int | float]]:
                     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]))
+                    # 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))
+                    # 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)
+                    lat,lon, _ = ctr.TransformPoint(e, n)
                     xyz.append([float(lat), float(lon), unc])
                 i += ir