|
| 1 | +# /// script |
| 2 | +# requires-python = ">=3.12" |
| 3 | +# dependencies = [ |
| 4 | +# "bedrock-ge==0.2.3", |
| 5 | +# "chardet==5.2.0", |
| 6 | +# "folium==0.19.5", |
| 7 | +# "geopandas==1.0.1", |
| 8 | +# "mapclassify==2.8.1", |
| 9 | +# "marimo", |
| 10 | +# "matplotlib==3.10.1", |
| 11 | +# "pandas==2.2.3", |
| 12 | +# "pyproj==3.7.1", |
| 13 | +# "requests==2.32.3", |
| 14 | +# "shapely==2.1.0", |
| 15 | +# "pygef"==0.11.1" |
| 16 | +# ] |
| 17 | +# /// |
| 18 | + |
| 19 | +import marimo |
| 20 | + |
| 21 | +__generated_with = "0.13.2" |
| 22 | +app = marimo.App(width="medium") |
| 23 | + |
| 24 | + |
| 25 | +@app.cell |
| 26 | +def _(mo): |
| 27 | + mo.md( |
| 28 | + """ |
| 29 | + # GEF Data for A2 Tunnel Maastricht |
| 30 | +
|
| 31 | + This notebook demonstrates how to |
| 32 | +
|
| 33 | + 1. Read in GEF files using [pygef](https://cemsbv.github.io/pygef/) |
| 34 | + 1. Use `bedrock-ge` to load Ground Investigation (GI) data from these GEF files |
| 35 | + 1. Convert that data into a standardized GI database using `bedrock-ge` |
| 36 | + 1. Transform the GI data into 3D GIS features with proper coordinates and geometry ([OGC Simple Feature Access](https://en.wikipedia.org/wiki/Simple_Features)) |
| 37 | + 1. Explore and analyze the GI data using: |
| 38 | + * Interactive filtering with Pandas dataframes |
| 39 | + * Visualization on interactive maps with GeoPandas |
| 40 | + 1. Export the processed GI database to a GeoPackage file for use in GIS software |
| 41 | +
|
| 42 | + ## Context |
| 43 | +
|
| 44 | + The Koning Willem-Alexander Tunnel is a double-deck tunnel for motorized traffic in the Dutch city of Maastricht. The tunnel has a length of 2.5 kilometers (lower tunnel tubes) and 2.3 kilometers (upper tunnel tubes). |
| 45 | +
|
| 46 | + The tunnel has moved the old A2 highway underground. It previously formed a barrier for the city and slowed traffic. |
| 47 | +
|
| 48 | + ### Geology |
| 49 | +
|
| 50 | +
|
| 51 | + The soil here consists of a coarse gravel layer with limestone limestone with karst phenomena underneath. There are also faults. |
| 52 | + Zwerfkeien |
| 53 | +
|
| 54 | + [Geotechniek-en-Risicos bij A2 Maastricht](https://www.cob.nl/magazines-brochures-en-nieuws/verdieping/verdieping-sept2012/geotechniek-en-risicos-bij-a2-maastricht/) |
| 55 | + [Source](https://archisarchief.cultureelerfgoed.nl/Archis2/Archeorapporten/24/AR26905/RAP%202709_4130060%20Maastricht%20A2-traverse.pdf) |
| 56 | +
|
| 57 | + ## Ground Investigation Data |
| 58 | +
|
| 59 | + The GI data was downloaded from [Dinoloket](https://www.dinoloket.nl/ondergrondgegevens), a platform where you can view and request data and models from TNO and BRO about the subsurface of the Netherlands. |
| 60 | + """ |
| 61 | + ) |
| 62 | + return |
| 63 | + |
| 64 | + |
| 65 | +@app.cell |
| 66 | +def _(Path, pygef): |
| 67 | + folder_path = Path("./gefs") |
| 68 | + gef_files = list(folder_path.glob("*.gef")) |
| 69 | + boreholes = [pygef.read_bore(gef_file) for gef_file in gef_files] |
| 70 | + return (boreholes,) |
| 71 | + |
| 72 | + |
| 73 | +@app.cell(hide_code=True) |
| 74 | +def _(mo): |
| 75 | + mo.md("pygef uses [polars](https://pola.rs/) for DataFrames, for consistency we will convert them to [Pandas](http://pandas.pydata.org/) DataFrames in this notebook.").callout("warn") |
| 76 | + return |
| 77 | + |
| 78 | + |
| 79 | +@app.cell(hide_code=True) |
| 80 | +def _(boreholes, mo): |
| 81 | + options = {d.alias: i for i, d in enumerate(boreholes)} |
| 82 | + multiselect = mo.ui.dropdown(options, label="Select borehole") |
| 83 | + multiselect |
| 84 | + return (multiselect,) |
| 85 | + |
| 86 | + |
| 87 | +@app.cell(hide_code=True) |
| 88 | +def _(multiselect): |
| 89 | + index = multiselect.value or 0 |
| 90 | + return (index,) |
| 91 | + |
| 92 | + |
| 93 | +@app.cell(hide_code=True) |
| 94 | +def _(boreholes, index): |
| 95 | + boreholes[index].data.to_pandas().dropna(axis=1, how='all') # drop empty columns for display |
| 96 | + return |
| 97 | + |
| 98 | + |
| 99 | +@app.cell |
| 100 | +def _(boreholes, index, pygef): |
| 101 | + pygef.plotting.plot_bore(boreholes[index]) |
| 102 | + return |
| 103 | + |
| 104 | + |
| 105 | +@app.cell(hide_code=True) |
| 106 | +def _(mo): |
| 107 | + mo.md( |
| 108 | + r""" |
| 109 | + ## Converting multiple GEF files to a relational database |
| 110 | +
|
| 111 | + First, let's check in which projected coordinate system the provided data was recorded |
| 112 | + """ |
| 113 | + ) |
| 114 | + return |
| 115 | + |
| 116 | + |
| 117 | +@app.cell |
| 118 | +def _(boreholes): |
| 119 | + boreholes[1].delivered_location.srs_name |
| 120 | + return |
| 121 | + |
| 122 | + |
| 123 | +@app.cell |
| 124 | +def _(boreholes): |
| 125 | + code = {bore.delivered_location.srs_name for bore in boreholes}.pop() |
| 126 | + # epsg_code = code.split("EPSG::")[-1] |
| 127 | + # crs = f"EPSG:{epsg_code}" |
| 128 | + crs = "EPSG:7415" |
| 129 | + return (crs,) |
| 130 | + |
| 131 | + |
| 132 | +@app.cell(hide_code=True) |
| 133 | +def _(mo): |
| 134 | + mo.md( |
| 135 | + r""" |
| 136 | + EPSG:28992 is [Rijksdriehoekscoördinaten](https://nl.wikipedia.org/wiki/Rijksdriehoeksco%C3%B6rdinaten), also called Amersfoort / RD New. |
| 137 | +
|
| 138 | + EPSG:7415 is with Amersfoort / RD New + NAP height. ie with elevation which we need for 3D geometry |
| 139 | + """ |
| 140 | + ) |
| 141 | + return |
| 142 | + |
| 143 | + |
| 144 | +@app.cell |
| 145 | +def _(): |
| 146 | + wgs = "EPSG:4326" |
| 147 | + return |
| 148 | + |
| 149 | + |
| 150 | +@app.cell |
| 151 | +def _(): |
| 152 | + project_uid = "Maastricht A2" |
| 153 | + return (project_uid,) |
| 154 | + |
| 155 | + |
| 156 | +@app.cell |
| 157 | +def _(CRS, crs, pd, project_uid): |
| 158 | + project = pd.DataFrame({ |
| 159 | + "project_uid": [project_uid], # primary key |
| 160 | + "crs_wkt": CRS(crs).to_wkt() |
| 161 | + }) |
| 162 | + return (project,) |
| 163 | + |
| 164 | + |
| 165 | +@app.cell |
| 166 | +def _(insitu_geo, locations, project): |
| 167 | + brgi_db = {"Project": project, "Location": locations.drop(columns=["data"]), "InSitu_GEOL": insitu_geo } |
| 168 | + return (brgi_db,) |
| 169 | + |
| 170 | + |
| 171 | +@app.function |
| 172 | +def process_data(bore): |
| 173 | + df = bore.data.to_pandas().dropna(axis=1, how='all').rename(columns= |
| 174 | + { |
| 175 | + 'upperBoundary': 'depth_to_top', |
| 176 | + 'lowerBoundary': 'depth_to_base', |
| 177 | + 'upperBoundaryOffset': 'elevation_at_top', |
| 178 | + 'lowerBoundaryOffset': 'elevation_at_base' |
| 179 | + }) |
| 180 | + |
| 181 | + return df |
| 182 | + |
| 183 | + |
| 184 | +@app.cell |
| 185 | +def _(boreholes, pd, project_uid): |
| 186 | + locations_df = pd.DataFrame([ |
| 187 | + { |
| 188 | + "location_uid": f"{borehole.alias} {project_uid}", # primary key |
| 189 | + "project_uid": project_uid, # foreign key |
| 190 | + "data": process_data(borehole), |
| 191 | + "location_source_id": borehole.alias, |
| 192 | + "date": borehole.research_report_date, |
| 193 | + "location_type": "Hole", |
| 194 | + "easting": borehole.delivered_location.x, |
| 195 | + "northing": borehole.delivered_location.y, |
| 196 | + "depth_to_base": min(borehole.data["lowerBoundaryOffset"]), |
| 197 | + "ground_level_elevation": borehole.delivered_vertical_position_offset, |
| 198 | + "elevation_at_base": borehole.delivered_vertical_position_offset - min(borehole.data["lowerBoundaryOffset"]), |
| 199 | + } |
| 200 | + for borehole in boreholes |
| 201 | + ]) |
| 202 | + return (locations_df,) |
| 203 | + |
| 204 | + |
| 205 | +@app.cell |
| 206 | +def _(calculate_location_gis_geometry, crs, locations_df): |
| 207 | + locations = calculate_location_gis_geometry(locations_df, crs=crs) |
| 208 | + return (locations,) |
| 209 | + |
| 210 | + |
| 211 | +@app.cell |
| 212 | +def _(create_lon_lat_height_table, crs, locations): |
| 213 | + create_lon_lat_height_table(locations, crs).explore() |
| 214 | + return |
| 215 | + |
| 216 | + |
| 217 | +@app.cell |
| 218 | +def _(locations, pd): |
| 219 | + insitu = pd.DataFrame([ |
| 220 | + { |
| 221 | + **layer, |
| 222 | + "location_uid": location["location_uid"], # foreignkey |
| 223 | + "project_uid": location["project_uid"], # foreignkey |
| 224 | + } |
| 225 | + # Outer loop: iterate through each location |
| 226 | + for location in locations.to_dict('records') |
| 227 | + # Inner loop: iterate through each layer in the location's data dataframe |
| 228 | + for layer in location["data"].to_dict('records') |
| 229 | + ]) |
| 230 | + return (insitu,) |
| 231 | + |
| 232 | + |
| 233 | +@app.cell |
| 234 | +def _(insitu): |
| 235 | + insitu |
| 236 | + return |
| 237 | + |
| 238 | + |
| 239 | +@app.cell |
| 240 | +def _(calculate_in_situ_gis_geometry, crs, insitu, locations): |
| 241 | + insitu_geo = calculate_in_situ_gis_geometry(insitu, locations, crs) |
| 242 | + return (insitu_geo,) |
| 243 | + |
| 244 | + |
| 245 | +@app.cell |
| 246 | +def _(insitu_geo): |
| 247 | + insitu_geo |
| 248 | + return |
| 249 | + |
| 250 | + |
| 251 | +@app.cell(hide_code=True) |
| 252 | +def _(mo): |
| 253 | + mo.md( |
| 254 | + r""" |
| 255 | + ## Saving the GI geospatial database as a GeoPackage (.gpkg) |
| 256 | +
|
| 257 | + Finally, lets write it to an actual geospatial database file, so we can share our GI data with others. For example, to reuse it in other notebooks, create dashboards, access the GI data in QGIS or ArcGIS, and more... |
| 258 | +
|
| 259 | + A GeoPackage is an <abbr title="Open Geospatial Consortium">OGC-standardized</abbr> extension of SQLite (a relational database in a single file, .sqlite or .db) that allows you to store any type of GIS data (both raster as well as vector data) in a single file that has the .gpkg extension. Therefore, many (open-source) GIS software packages support GeoPackage! |
| 260 | + """ |
| 261 | + ) |
| 262 | + return |
| 263 | + |
| 264 | + |
| 265 | +@app.cell |
| 266 | +def _(brgi_db, check_brgi_database): |
| 267 | + check_brgi_database(brgi_db) |
| 268 | + return |
| 269 | + |
| 270 | + |
| 271 | +@app.cell |
| 272 | +def _(brgi_db, write_gi_db_to_gpkg): |
| 273 | + write_gi_db_to_gpkg(brgi_db, gpkg_path="./output/A2_Maastricht.gpkg") |
| 274 | + return |
| 275 | + |
| 276 | + |
| 277 | +@app.cell |
| 278 | +def _(): |
| 279 | + import marimo as mo |
| 280 | + import pygef |
| 281 | + import os |
| 282 | + from pathlib import Path |
| 283 | + import pandas as pd |
| 284 | + import geopandas as gpd |
| 285 | + import matplotlib |
| 286 | + import pyarrow |
| 287 | + import folium |
| 288 | + import mapclassify |
| 289 | + from shapely.geometry import Point, LineString |
| 290 | + from bedrock_ge.gi.gis_geometry import calculate_wgs84_coordinates, calculate_in_situ_gis_geometry, calculate_gis_geometry, calculate_location_gis_geometry, create_lon_lat_height_table |
| 291 | + from bedrock_ge.gi.write import write_gi_db_to_gpkg |
| 292 | + from bedrock_ge.gi.validate import check_brgi_database |
| 293 | + from pyproj import CRS |
| 294 | + return ( |
| 295 | + CRS, |
| 296 | + Path, |
| 297 | + calculate_in_situ_gis_geometry, |
| 298 | + calculate_location_gis_geometry, |
| 299 | + check_brgi_database, |
| 300 | + create_lon_lat_height_table, |
| 301 | + mo, |
| 302 | + pd, |
| 303 | + pygef, |
| 304 | + write_gi_db_to_gpkg, |
| 305 | + ) |
| 306 | + |
| 307 | + |
| 308 | +if __name__ == "__main__": |
| 309 | + app.run() |
0 commit comments