Create an API to allow any input GI to be mapped to the Bedrock GI Database model

.github/workflows/python-tests.yml

@@ -12,7 +12,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, windows-latest]
-        python-version: ['3.9', '3.12']
+        python-version: ['3.10', '3.12']
 
     steps:
     - name: Checkout code

.gitignore

@@ -6,6 +6,8 @@ __marimo__/
 **/log
 **/tmp
 **/temp
+examples/**/raw/
+examples/**/cleaned/
 
 # SQLite Databases
 *.sqlite

.python-version

@@ -0,0 +1 @@
+3.10

pyproject.toml

@@ -7,7 +7,7 @@ authors = [
 ]
 license = {text = "Apache Software License (Apache 2.0)"}
 readme = "README.md"
-requires-python = ">=3.9"
+requires-python = ">=3.10"
 dependencies = [
     "geopandas~=1.0",
     "openpyxl~=3.0",
@@ -55,7 +55,6 @@ classifiers = [
     "Topic :: Scientific/Engineering",
     "Topic :: Scientific/Engineering :: GIS",
     "Programming Language :: Python",
-    "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",

sandbox/pyproj_3d_transformations.py

@@ -0,0 +1,100 @@
+import marimo
+
+__generated_with = "0.13.11"
+app = marimo.App(width="medium")
+
+
+@app.cell
+def _():
+    import marimo as mo
+    from pyproj import CRS, Transformer
+    from pyproj.crs.crs import CompoundCRS
+    return CRS, CompoundCRS, Transformer
+
+
+@app.cell
+def _(CRS, Transformer):
+    # https://pyproj4.github.io/pyproj/stable/advanced_examples.html#promote-crs-to-3d
+    wgs84 = CRS("EPSG:4326")
+    swiss_proj = CRS("EPSG:2056")
+    transformer = Transformer.from_crs(wgs84, swiss_proj, always_xy=True)
+    # 2D Transformation
+    print(f"2D transform of point {(8.37909, 47.01987, 1000)} from {wgs84} to {swiss_proj} gives:")
+    print(transformer.transform(8.37909, 47.01987, 1000))
+
+    wgs84_3d = wgs84.to_3d()
+    swiss_proj_ellipsoidal_3d = swiss_proj.to_3d()
+    transformer_ellipsoidal_3d = Transformer.from_crs(
+        wgs84_3d,
+        swiss_proj_ellipsoidal_3d,
+        always_xy=True,
+    )
+    # 3D Transformation
+    print(f"3D transform of point {(8.37909, 47.01987, 1000)} from {wgs84_3d.to_string()} to {swiss_proj_ellipsoidal_3d} gives:")
+    print(transformer_ellipsoidal_3d.transform(8.37909, 47.01987, 1000))
+
+    return swiss_proj, transformer_ellipsoidal_3d, wgs84_3d
+
+
+@app.cell
+def _(
+    CRS,
+    CompoundCRS,
+    Transformer,
+    swiss_proj,
+    transformer_ellipsoidal_3d,
+    wgs84_3d,
+):
+    # https://pyproj4.github.io/pyproj/stable/build_crs.html#compound-crs
+    swiss_lhn95_height = CRS("EPSG:5729")
+    swiss_compound = CompoundCRS(
+        name="CH1903+ / LV95 + LHN95 height",
+        components=[swiss_proj, swiss_lhn95_height]
+    )
+    transformer_wgs84_3d_to_swiss_compound = Transformer.from_crs(
+        wgs84_3d,
+        swiss_compound,
+        always_xy=True,
+    )
+    print(transformer_ellipsoidal_3d.transform(8.37909, 47.01987, 1000))
+    return
+
+
+@app.cell
+def _(CRS):
+    uk_grid_27700 = CRS(27700)
+    uk_3d = uk_grid_27700.to_3d()
+    swiss_2056 = CRS("EPSG:2056")
+    swiss_3d = swiss_2056.to_3d()
+    egm2008_3855 = CRS(3855)
+    rdnew_nap_7415 = CRS(7415)
+    wgs84_egm2008_9518 = CRS(9518)
+    crs_components = extract_crs_components(swiss_3d)
+    crs_components
+    return
+
+
+@app.function
+def extract_crs_components(compound_crs):
+    """Extract horizontal and vertical CRS from a compound CRS"""
+    if not compound_crs.is_compound:
+        return compound_crs, None
+
+    horizontal_crs = None
+    vertical_crs = None
+
+    for sub_crs in compound_crs.sub_crs_list:
+        if sub_crs.is_projected or sub_crs.is_geographic:
+            print(f"Horizontal CRS {sub_crs.name} is a {sub_crs.type_name} and has EPSG:{sub_crs.to_epsg()}.")
+            horizontal_crs = sub_crs
+        elif sub_crs.is_vertical:
+            print(f"Vertical CRS {sub_crs.name} has EPSG:{sub_crs.to_epsg()}.")
+            vertical_crs = sub_crs
+        else:
+            print(f"This CRS is not horizontal (projected or geographic) nor vertical: {sub_crs.type_name}")
+
+    return horizontal_crs, vertical_crs
+
+
+if __name__ == "__main__":
+    app.run()

src/bedrock_ge/gi/ags.py

@@ -0,0 +1,103 @@
+from pathlib import Path
+from typing import IO
+
+from pyproj import CRS
+
+from bedrock_ge.gi.ags3 import ags3_to_brgi_db_mapping
+from bedrock_ge.gi.io_utils import detect_encoding, open_text_data_source
+from bedrock_ge.gi.mapping_models import BedrockGIMapping
+
+
+def ags_to_brgi_db_mapping(
+    source: str | Path | IO[str] | IO[bytes] | bytes,
+    projected_crs: CRS,
+    vertical_crs: CRS = CRS(3855),
+    encoding: str | None = None,
+) -> BedrockGIMapping:
+    """Map AGS 3 or AGS 4 data to the Bedrock GI data model.
+
+    Args:
+        source (str | Path | IO[str] | IO[bytes] | bytes): The AGS file (str or Path)
+            or a file-like object that represents the AGS file.
+        projected_crs (CRS): Projected Coordinate Reference System (CRS). For example:
+            - OSGB36 / British National Grid: `pyproj.CRS("EPSG:27700")`
+            - Hong Kong 1980 Grid System: `pyproj.CRS("EPSG:2326")`
+        vertical_crs (CRS, optional): Vertical CRS. Defaults to EGM2008 height, EPSG:3855
+            which measures the orthometric height w.r.t. the Earth Gravitational Model 2008.
+            - Ordnance Datum Newlyn (ODN) Height: `pyproj.CRS("EPSG:5701")`
+            - Hong Kong Principle Datum (HKPD) Height: `pyproj.CRS("EPSG:5738")`
+        encoding (str | None, optional): Encoding of the text file or bytes stream.
+            Defaults to None. An attempt at detecting the encoding will be made if None.
+
+    Raises:
+        ValueError: If the data does not match AGS 3 or AGS 4 format.
+
+    Returns:
+        BedrockGIDatabaseMapping: Object that maps AGS 3 or AGS 4 data to Bedrock GI data model.
+    """
+    if not encoding:
+        encoding = detect_encoding(source)
+
+    # Get first non-blank line, None if all lines are blank
+    with open_text_data_source(source, encoding=encoding) as f:
+        first_line = next((line.strip() for line in f if line.strip()), None)
+
+    if first_line:
+        if first_line.startswith('"**'):
+            ags_version = 3
+            brgi_db_mapping = ags3_to_brgi_db_mapping(
+                source, projected_crs, vertical_crs, encoding
+            )
+        elif first_line.startswith('"GROUP"'):
+            ags_version = 4
+            # brgi_db_mapping = ags4_to_brgi_db_mapping(
+            #     source, projected_crs, vertical_crs, encoding
+            # )
+        else:
+            # If first non-empty line doesn't match AGS 3 or AGS 4 format
+            raise ValueError("The data provided is not valid AGS 3 or AGS 4 data.")
+    else:
+        raise ValueError("The file provided has only blank lines")
+
+    # Put CPT data into brgi_db.Other table, because CPT data has too many rows
+    # that would generate geospatial geometry.
+    # "STCN" and "SCPT" are the group names for CPT data in AGS 3 and AGS 3 respectively.
+    # TODO: implement a warning when interpolating GI geospatial geometry when
+    # TODO: a single GI location has waaay too many rows in a certain In-Situ test,
+    # TODO: rather than removing specific groups here.
+    insitu_test_names = {
+        insitu_test.table_name: i
+        for i, insitu_test in enumerate(brgi_db_mapping.InSitu)
+    }
+    cpt_key = (
+        "STCN"
+        if "STCN" in insitu_test_names
+        else "SCPT"
+        if "SCPT" in insitu_test_names
+        else None
+    )
+    if cpt_key is not None:
+        cpt_data_mapping = brgi_db_mapping.InSitu.pop(insitu_test_names[cpt_key])
+        del insitu_test_names[cpt_key]
+        brgi_db_mapping.Other.append(cpt_data_mapping)
+
+    # Log information about the mapped AGS 3 or AGS 4 data
+    project_id = brgi_db_mapping.Project.project_id
+    n_gi_locations = len(brgi_db_mapping.Location.data)
+    n_samples = len(brgi_db_mapping.Sample.data) if brgi_db_mapping.Sample else 0
+    print_args = [
+        f"AGS {ags_version} data was read for Project {project_id}",
+        f"This GI data contains {n_gi_locations} GI locations, {n_samples} samples and:",
+        f"  - In-Situ Tests: {list(insitu_test_names.keys())}",
+    ]
+    if brgi_db_mapping.Lab:
+        print_args.append(
+            f"  - Lab Tests: {[lab_test.table_name for lab_test in brgi_db_mapping.Lab]}"
+        )
+    if brgi_db_mapping.Other:
+        print_args.append(
+            f"  - Other Tables: {[other_table.table_name for other_table in brgi_db_mapping.Other]}"
+        )
+    print(*print_args, sep="\n", end="\n\n")
+
+    return brgi_db_mapping