Rename marimo notebook such that it doesn't conflict with Jules' CPT interpretation notebook

Author: JoostGevaert

examples/nl_amsterdam/nl_amsterdam_bro.py …amples/nl_amsterdam/get_cpts_from_bro.pyexamples/nl_amsterdam/nl_amsterdam_bro.py renamed to examples/nl_amsterdam/get_cpts_from_bro.py examples/nl_amsterdam/nl_amsterdam_bro.py renamed to examples/nl_amsterdam/get_cpts_from_bro.py

@@ -13,10 +13,167 @@
 import marimo
 
 __generated_with = "0.17.0"
-app = marimo.App(width="medium")
+app = marimo.App(width="columns")
+
+
+@app.cell(column=0, hide_code=True)
+def _(mo):
+    mo.md(
+        r"""
+    How to access BRO data: [Handreiking Afname BRO Gegevens](https://www.bro-productomgeving.nl/bpo/latest/handreiking-afname-bro-gegevens)
+
+    1. [BROloket](https://www.broloket.nl/ondergrondgegevens)
+    2. BRO APIs
+      1. SOAP - impractical, because requires a digital "PKI" certificate.
+      2. [REST](https://www.bro-productomgeving.nl/bpo/latest/url-s-publieke-rest-services)
+        1. CPT: <https://publiek.broservices.nl/sr/cpt/v1>
+        2. Goetechnical boreholes: <https://publiek.broservices.nl/sr/bhrgt/v2>
+    3. [PDOK](https://app.pdok.nl/viewer)
+      1. WMS - This is useful for quickly viewing the location of historic CPTs or geotechnical boreholes.
+      2. ATOM feed - For downloading the whole dataset, i.e. to download all CPTs or all geotechnical boreholes in BRO.
+    """
+    )
+    return
+
+
+@app.cell(hide_code=True)
+def _(Request, bounds, datetime, json, minidom, mo, pl, urlopen, xmltodict):
+    cpt_search_url = "https://publiek.broservices.nl/sr/cpt/v1/characteristics/searches"
+    api_request_data = json.dumps(
+        {
+            "registrationPeriod": {
+                "beginDate": "2017-01-01",
+                "endDate": datetime.date.today().isoformat(),
+            },
+            "area": {
+                "boundingBox": {
+                    "lowerCorner": {
+                        "lon": bounds[0],
+                        "lat": bounds[1],
+                    },
+                    "upperCorner": {
+                        "lon": bounds[2],
+                        "lat": bounds[3],
+                    },
+                }
+            },
+        }
+    ).encode("utf-8")
+
+    cpt_search_req = Request(
+        cpt_search_url,
+        data=api_request_data,
+        headers={"Content-Type": "application/json"},
+        method="POST",
+    )
+    with urlopen(cpt_search_req, timeout=30) as cpt_search_resp:
+        xml = cpt_search_resp.read()
+
+    cpt_search_dict = xmltodict.parse(
+        xml,
+        xml_attribs=False,
+        process_namespaces=True,
+        namespaces={
+            "http://www.opengis.net/gml/3.2": None,
+            "http://www.broservices.nl/xsd/dscpt/1.1": None,
+            "http://www.broservices.nl/xsd/brocommon/3.0": None,
+        },
+    )
+
+    def to_bool(col: str) -> pl.Expr:
+        truthy = ["ja", "yes", "true", "1"]
+        return pl.col(col).str.to_lowercase().is_in(truthy)
+
+    def parse_pos(col: str) -> pl.Expr:
+        return (
+            pl.col(col)
+            .struct.field("pos")
+            .str.split(" ")
+            .cast(pl.Array(float, 2))
+            .alias(col)
+        )
+
+    cpt_search_df = (
+        pl.DataFrame(
+            cpt_search_dict["dispatchCharacteristicsResponse"]["dispatchDocument"]
+        )
+        .unnest("CPT_C")
+        .with_columns(
+            # Casts
+            pl.col("deliveryAccountableParty").cast(int),
+            pl.col("objectRegistrationTime").cast(pl.Datetime),
+            pl.col("offset").cast(float),
+            pl.col("startTime").cast(pl.Date),
+            pl.col("predrilledDepth").cast(float),
+            pl.col("finalDepth").cast(float),
+            # "ja" / "nee" -> boolean
+            to_bool("deregistered"),
+            to_bool("underReview"),
+            to_bool("dissipationTestPerformed"),
+            # {"pos":"lon lat"} -> pl.Array[float,2 ] (keeps [lon, lat])
+            parse_pos("standardizedLocation"),
+            parse_pos("deliveredLocation"),
+            # {"date":"YYYY-MM-DD"} -> pl.Date
+            pl.col("researchReportDate")
+            .struct.field("date")
+            .cast(pl.Date)
+            .alias("researchReportDate"),
+        )
+    )
+    cpt_search_table = mo.ui.table(cpt_search_df)
+
+    tabs = mo.ui.tabs(
+        {
+            "DataFrame": cpt_search_table,
+            "XML": mo.md(
+                f"```xml\n{minidom.parseString(xml).toprettyxml(indent='  ')}\n```"
+            ),
+            "JSON": xmltodict.parse(xml),
+            "JSON, no XML attributes": cpt_search_dict,
+        }
+    )
+
+    tabs
+    return (cpt_search_df,)
 
 
 @app.cell
+def _(BytesIO, Request, cpt_search_df, pl, read_cpt, urlopen):
+    cpts = []
+    cpt_data = []
+    for bro_id in cpt_search_df["broId"]:
+        cpt_get_url = f"https://publiek.broservices.nl/sr/cpt/v1/objects/{bro_id}"
+        cpt_get_req = Request(cpt_get_url, method="GET")
+        with urlopen(cpt_get_req, timeout=30) as resp:
+            cpt_xml = resp.read()
+            cpt_obj = read_cpt(BytesIO(cpt_xml))
+        cpt_data.append(
+            cpt_obj.data.with_columns(pl.lit(cpt_obj.bro_id).alias("broId"))
+        )
+        cpt_dict = cpt_obj.__dict__.copy()
+        del cpt_dict["data"]
+        cpt_dict["delivered_location"] = [
+            cpt_dict["delivered_location"].x,
+            cpt_dict["delivered_location"].y,
+        ]
+        cpt_dict["standardized_location"] = [
+            cpt_dict["standardized_location"].x,
+            cpt_dict["standardized_location"].y,
+        ]
+        cpt_dict["delivered_vertical_position_datum"] = cpt_dict["delivered_vertical_position_datum"]._name_
+        cpts.append(cpt_dict)
+
+    cpts_df = pl.DataFrame(cpts).cast(
+        {
+            "delivered_location": pl.Array(float, 2),
+            "standardized_location": pl.Array(float, 2),
+        }
+    )
+    cpt_data_df = pl.concat(cpt_data, how="diagonal_relaxed")
+    return (cpt_data_df,)
+
+
+@app.cell(hide_code=True)
 def _():
     import datetime
     import json
@@ -27,7 +184,9 @@ def _():
     import folium
     import geopandas as gpd
     import marimo as mo
+    import polars as pl
     import xmltodict
+
     from folium.plugins import Draw
     from lxml import etree
     from pygef import read_cpt
@@ -39,40 +198,19 @@ def _():
         Draw,
         Request,
         datetime,
-        etree,
         folium,
         gpd,
         json,
         minidom,
         mo,
-        plot_cpt,
+        pl,
         read_cpt,
         urlopen,
         xmltodict,
     )
 
 
-@app.cell
-def _(mo):
-    mo.md(
-        r"""
-    How to access BRO data: [Handreiking Afname BRO Gegevens](https://www.bro-productomgeving.nl/bpo/latest/handreiking-afname-bro-gegevens)
-
-    1. [BROloket](https://www.broloket.nl/ondergrondgegevens)
-    2. BRO APIs
-      1. SOAP - impractical, because requires a digital "PKI" certificate.
-      2. [REST](https://www.bro-productomgeving.nl/bpo/latest/url-s-publieke-rest-services)
-        1. CPT: <https://publiek.broservices.nl/sr/cpt/v1>
-        2. Goetechnical boreholes: <https://publiek.broservices.nl/sr/bhrgt/v2>
-    3. [PDOK](https://app.pdok.nl/viewer)
-      1. WMS - This is useful for quickly viewing the location of historic CPTs or geotechnical boreholes.
-      2. ATOM feed - For downloading the whole dataset, i.e. to download all CPTs or all geotechnical boreholes in BRO.
-    """
-    )
-    return
-
-
-@app.cell(hide_code=True)
+@app.cell(column=1, hide_code=True)
 def _(Draw, buffer, folium, geojson_text_area, gpd, json, site_geojson):
     geojson = geojson_text_area.value if geojson_text_area.value else site_geojson
     # Create a folium interactive map (leaflet.js maps)
@@ -87,16 +225,18 @@ def _(Draw, buffer, folium, geojson_text_area, gpd, json, site_geojson):
     buffered_site = site.geometry.buffer(buffer.value)
     bounds = buffered_site.to_crs(4326).bounds.to_numpy()[0]
     folium_map = buffered_site.explore(
-        tiles="https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}",
+        # tiles="https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}",
+        tiles="CartoDB positron",
         style_kwds={"fillOpacity": 0.1},
         attr=("Esri.WorldStreetMap"),
     )
     site.explore(m=folium_map, color="red", style_kwds={"fill": False})
-    folium.Rectangle(
+    bounds_rectangle = folium.Rectangle(
         bounds=[[bounds[1], bounds[0]], [bounds[3], bounds[2]]],
         color="black",
         weight=0.5,
-    ).add_to(folium_map)
+    )
+    bounds_rectangle.add_to(folium_map)
 
     # Add PDOK's CPT WMS layer
     folium.WmsTileLayer(
@@ -132,10 +272,10 @@ def _(Draw, buffer, folium, geojson_text_area, gpd, json, site_geojson):
     ).add_to(folium_map)
 
     folium_map
-    return (bounds,)
+    return bounds, bounds_rectangle, site
 
 
-@app.cell
+@app.cell(hide_code=True)
 def _(mo):
     site_geojson = '{"type":"Feature","properties":{},"geometry":{"type":"Polygon","coordinates":[[[4.902889698063478,52.385435257999404],[4.902923632305726,52.385509906165986],[4.902901716440955,52.385619936470341],[4.903006347021203,52.38569026942325],[4.90389499998988,52.385854235381281],[4.904091535809526,52.385902562073795],[4.904285950739031,52.385988859607401],[4.904484607448818,52.38602683046868],[4.904532680958658,52.386042363993511],[4.904790722592364,52.386094573856099],[4.904799206152918,52.386073862513044],[4.904974533071164,52.386108812898868],[4.904999276789463,52.386100183176488],[4.905287541107636,52.385853480276289],[4.903576336246269,52.385074312916224],[4.903554420381488,52.385093082996256],[4.903498923756159,52.38509416173622],[4.903484784488561,52.385100202679702],[4.903473119592793,52.385110990076754],[4.903287895187223,52.385027927051986],[4.903256788798504,52.38504302943177],[4.903150390809821,52.38513472234078],[4.903119637902792,52.385127171167227],[4.903098782483085,52.38514421524286],[4.903084996697174,52.385139037296511],[4.90302349088312,52.385198152147602],[4.902889698063478,52.385435257999404]]]}}'
     geojson_text_area = mo.ui.text_area(
@@ -151,108 +291,37 @@ def _(mo):
     return buffer, geojson_text_area, site_geojson
 
 
-@app.cell
-def _(buffer):
-    # EPSG:4258 is the European Terrestrial Reference System 1989 (ETR)
-    # ETR89 coordinates are in [Longitude (x), Latitude (y)]
-    lonlat_bounds = buffer.to_crs(4258).bounds
-    lonlat_bounds
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""The bouding box around the site + buffer that we created above is now used to make a call to the BRO API to retrieve all the CPT characteristics of the CPTs in that bounding box.""")
     return
 
 
-@app.cell
-def _(Request, bounds, datetime, json, minidom, mo, urlopen, xmltodict):
-    cpt_search_url = "https://publiek.broservices.nl/sr/cpt/v1/characteristics/searches"
-    api_request_data = json.dumps(
-        {
-            "registrationPeriod": {
-                "beginDate": "2017-01-01",
-                "endDate": datetime.date.today().isoformat(),
-            },
-            "area": {
-                "boundingBox": {
-                    "lowerCorner": {
-                        "lon": bounds[0],
-                        "lat": bounds[1],
-                    },
-                    "upperCorner": {
-                        "lon": bounds[2],
-                        "lat": bounds[3],
-                    },
-                }
-            },
-        }
-    ).encode("utf-8")
-
-    cpt_search_req = Request(
-        cpt_search_url,
-        data=api_request_data,
-        headers={"Content-Type": "application/json"},
-        method="POST",
-    )
-    with urlopen(cpt_search_req, timeout=30) as cpt_search_resp:
-        xml = cpt_search_resp.read()
-
-    cpt_search = xmltodict.parse(
-        xml,
-        xml_attribs=False,
-        process_namespaces=True,
-        namespaces={
-            "http://www.broservices.nl/xsd/dscpt/1.1": None,
-            "http://www.broservices.nl/xsd/brocommon/3.0": None,
-        },
-    )
-
-    tabs = mo.ui.tabs(
-        {
-            "XML": mo.md(
-                f"```xml\n{minidom.parseString(xml).toprettyxml(indent='  ')}\n```"
-            ),
-            "JSON": xmltodict.parse(xml),
-            "JSON, no XML attributes": cpt_search,
-        }
-    )
-
-    tabs
+@app.cell(hide_code=True)
+def _(bounds_rectangle, cpt_search_df, gpd, site):
+    cpt_search_gdf = gpd.GeoDataFrame(cpt_search_df.to_pandas(), geometry=gpd.points_from_xy(cpt_search_df["deliveredLocation"].arr.get(0), cpt_search_df["deliveredLocation"].arr.get(1)), crs=28992)
+    cpt_map = cpt_search_gdf.explore(tiles="CartoDB positron")
+    bounds_rectangle.add_to(cpt_map)
+    site.explore(m=cpt_map, color="red", style_kwds={"fill": False})
+    cpt_map
     return
 
 
-app._unparsable_cell(
-    r"""
-    for cpt in cpt_search
-    """,
-    name="_",
-)
-
-
 @app.cell
-def _(Request, etree, urlopen):
-    cpt_get_url = "https://publiek.broservices.nl/sr/cpt/v1/objects/CPT000000198164"
-    cpt_get_req = Request(cpt_get_url, method="GET")
-    with urlopen(cpt_get_req, timeout=30) as resp:
-        cpt_str = resp.read()
-        cpt_xml = etree.fromstring(cpt_str)
-
-    print(etree.tostring(cpt_xml, pretty_print=True, encoding="unicode"))
-    return (cpt_str,)
-
-
-@app.cell
-def _(cpt_str):
-    cpt_str
+def _(cpt_data_df):
+    cpt_data_df
     return
 
 
 @app.cell
-def _(BytesIO, cpt_str, read_cpt):
-    cpt_data = read_cpt(BytesIO(cpt_str))
-    cpt_data.__dict__
-    return (cpt_data,)
+def _(cpt_data_df, mo):
+    mo.ui.dataframe(cpt_data_df)
+    return
 
 
 @app.cell
-def _(cpt_data, plot_cpt):
-    plot_cpt(cpt_data)
+def _(cpt_data_df):
+    cpt_data_df.filter(cpt_data_df["broId"].contains("CPT000000198163"))
     return