From 8d7c90682e9a341bba0d0ed067550092b78b9c85 Mon Sep 17 00:00:00 2001
From: Trey Stafford <trey.stafford@colorado.edu>
Date: Mon, 26 Jun 2023 15:43:56 -0600
Subject: [PATCH] Experiment with creating contours from GEBCO grid

---
 .../Bathymetry/bathymetric_chart.py           | 75 +++++++++++++++----
 1 file changed, 59 insertions(+), 16 deletions(-)

diff --git a/qgreenland/config/layers/Oceanography/Bathymetry/bathymetric_chart.py b/qgreenland/config/layers/Oceanography/Bathymetry/bathymetric_chart.py
index 0c47a19d..4ddc20d8 100644
--- a/qgreenland/config/layers/Oceanography/Bathymetry/bathymetric_chart.py
+++ b/qgreenland/config/layers/Oceanography/Bathymetry/bathymetric_chart.py
@@ -1,5 +1,4 @@
 from qgreenland.config.datasets.bathymetric_chart import gebco_bathymetric_chart
-from qgreenland.config.helpers.layers.geological_map import make_layer
 from qgreenland.config.helpers.steps.compress_and_add_overviews import (
     compress_and_add_overviews,
 )
@@ -7,6 +6,7 @@
 from qgreenland.config.helpers.steps.warp_and_cut import warp_and_cut
 from qgreenland.config.project import project
 from qgreenland.models.config.layer import Layer, LayerInput
+from qgreenland.models.config.step import CommandStep
 
 _background_boundary = project.boundaries["background"]
 _background_bbox = _background_boundary.bbox
@@ -16,6 +16,9 @@
     description=(
         """Bathymetric elevation in meters.
 
+        Negative values represent bathymetric depths and positive values are
+        topographic heights.
+
         Elevations can be assumed to be relative to mean sea level. However, in
         some shallow water areas, the grids include data from sources having a
         vertical datum other than mean sea level. We are working to understand
@@ -49,21 +52,61 @@
     ],
 )
 
+gebco_bathymetric_contours = Layer(
+    id="gebco_bathymetric_contours",
+    title="Depth contours",
+    description=(
+        """Bathymetric elevation contours in meters.
 
-bathymetric_contours_params = {
-    "id": "bathymetric_contours",
-    "title": "Depth contours",
-    "description": (
-        """This dataset includes linear features that represent
-         bathymetric contours recorded in metres,
-        derived from the International Bathymetric Chart of the Arctic Ocean."""
-    ),
-    "style": "bathymetry",
-    "input_filepath": "data/shape/base/bathymetry",
-    "fn_mask": "bathymetry.*",
-}
+        Contour interval is 100m and was derived from the "Depth (400m)" layer.
 
-bathymetric_contours = make_layer(
-    layer_id=bathymetric_contours_params["id"],
-    layer_params=bathymetric_contours_params,
+        Elevations can be assumed to be relative to mean sea level. However, in
+        some shallow water areas, the grids include data from sources having a
+        vertical datum other than mean sea level. We are working to understand
+        how best to fully assimilate these data."""
+    ),
+    tags=[],
+    style="bathymetry",
+    input=LayerInput(
+        dataset=gebco_bathymetric_chart,
+        asset=gebco_bathymetric_chart.assets["only"],
+    ),
+    steps=[
+        decompress_step(input_file="{input_dir}/gebco_2023.zip"),
+        *warp_and_cut(
+            input_file="NETCDF:{input_dir}/GEBCO_2023.nc:elevation",
+            output_file="{output_dir}/bathymetric_chart.tif",
+            reproject_args=(
+                "-te",
+                f"{_background_bbox.min_x} {_background_bbox.min_y} {_background_bbox.max_x} {_background_bbox.max_y}",
+                "-tr",
+                "400",
+                "400",
+            ),
+            cut_file=project.boundaries["background"].filepath,
+        ),
+        CommandStep(
+            args=[
+                "gdal_calc.py",
+                '--calc="where(A < 0, absolute(A), nan)"',
+                "--A_band=1",
+                "-A",
+                "{input_dir}/bathymetric_chart.tif",
+                "--outfile={output_dir}/positive_depths.tif",
+            ],
+        ),
+        CommandStep(
+            args=[
+                "gdal_contour",
+                "-a",
+                "DEPTH",
+                "-i",
+                "500",
+                "-nln",
+                "bathymetry",
+                "{input_dir}/positive_depths.tif",
+                "{output_dir}/bathymetric_chart.gpkg",
+            ],
+        ),
+    ],
 )