Implement location_gis_geometry function

Author: JoostGevaert

examples/hk_kaitak_ags3/hk_kaitak_ags3_to_brgi_geodb.py

@@ -291,10 +291,17 @@ def _():
 
 @app.cell
 def _(brgi_db):
-    from bedrock_ge.gi.geospatial import create_lon_lat_height_table
+    from bedrock_ge.gi.geospatial import create_lon_lat_height_table, location_gis_geometry
 
     lon_lat_height_gdf = create_lon_lat_height_table(brgi_db)
     lon_lat_height_gdf.explore()
+    return (location_gis_geometry,)
+
+
+@app.cell
+def _(brgi_db, location_gis_geometry):
+    location_gdf = location_gis_geometry(brgi_db)
+    location_gdf.explore()
     return
 
 

src/bedrock_ge/gi/geospatial.py

@@ -2,13 +2,50 @@
 import pandas as pd
 from pyproj import CRS, Transformer
 from pyproj.crs import CompoundCRS
+from shapely.geometry import LineString
 
 from bedrock_ge.gi.schemas import BedrockGIDatabase, BedrockGIGeospatialDatabase
 
 
 def location_gis_geometry(brgi_db: BedrockGIDatabase) -> gpd.GeoDataFrame:
-    location_gdf = brgi_db.Location
-    return location_gdf
+    # TODO: Implement logic to handle multiple CRS'es in the input GI data:
+    #       1. Create WKT geometry for each location in original CRS
+    #       2. Convert to WGS84 + EGM2008 orthometric height EPSG:9518
+    #       3. Interpolate InSituTest and Sample geospatial vector geometry from active geometry column
+    hor_crs = brgi_db.Project["horizontal_crs_wkt"]
+    vert_crs = brgi_db.Project["vertical_crs_wkt"]
+    if hor_crs.nunique() > 1 or vert_crs.nunique() > 1:
+        raise ValueError(
+            "All projects must have the same horizontal and vertical CRS (Coordinate Reference System).\n"
+            "Raise an issue on GitHub in case you need to be able to combine GI data that was acquired in multiple different CRSes."
+        )
+
+    hor_crs = CRS.from_wkt(hor_crs.iat[0])
+    vert_crs = CRS.from_wkt(vert_crs.iat[0])
+    compound_crs = CompoundCRS(
+        name=f"{hor_crs.name} + {vert_crs.name}",
+        components=[hor_crs, vert_crs],
+    )
+
+    # TODO: Implement logic such that inclined borholes are handled correctly.
+    #       All boreholes are now assumed to be vertical.
+    location_df = brgi_db.Location.copy()
+    location_df["elevation_at_base"] = (
+        location_df["ground_level_elevation"] - location_df["depth_to_base"]
+    )
+    return gpd.GeoDataFrame(
+        brgi_db.Location.copy(),
+        geometry=location_df.apply(
+            lambda row: LineString(
+                [
+                    (row["easting"], row["northing"], row["ground_level_elevation"]),
+                    (row["easting"], row["northing"], row["elevation_at_base"]),
+                ]
+            ),
+            axis=1,
+        ),
+        crs=compound_crs,
+    )
 
 
 def create_lon_lat_height_table(brgi_db: BedrockGIDatabase) -> gpd.GeoDataFrame:
@@ -45,10 +82,10 @@ def create_lon_lat_height_table(brgi_db: BedrockGIDatabase) -> gpd.GeoDataFrame:
     lon_lat_height_df = pd.concat(dfs, ignore_index=True)
     return gpd.GeoDataFrame(
         lon_lat_height_df,
-        crs=wgs84_egm2008_crs,
         geometry=gpd.points_from_xy(
             lon_lat_height_df["longitude"],
             lon_lat_height_df["latitude"],
             lon_lat_height_df["egm2008_ground_level_height"],
         ),
+        crs=wgs84_egm2008_crs,
     )