Trying with official WFSlib

01_baden_wuerttemberg.py

@@ -13,29 +13,97 @@ BADEN_WUERTTEMBERG = Bundesland(
 )
 
 if __name__ == "__main__":
-    bw_json = TMP_DIR / "bw_lod2.json"
-    download_file(
-        BW_WFS,
-        bw_json
-    )
-
-    with open(bw_json) as json_file:
-        data = json.load(json_file)
-    xs, ys = [], []
-    import matplotlib.pyplot as plt
-    for feature in data["features"]:
-        x, y = feature['bbox'][:2]
-        date = feature['properties']['produktionsdatum'] or '1970-01-01'
-        xs.append(x)
-        ys.append(y)
-        # citygmls.append(
-        #     CityGMLWithDate(
-        #         url=feature["properties"]["xml"],
-        #         bundesland=NIEDERSACHSEN,
-        #         source_date=datetime.strptime(feature["properties"]["Aktualitaet"], "%Y-%m-%d %H:%M:%S"),
-        #         # coordinate_reference_system ?
-        #     )
-        # )
-    # TODO: Show dates too?
-    plt.scatter(xs, ys, s=0.1)
-    plt.show()
+    # bw_json = TMP_DIR / "bw_lod2.json"
+    # download_file(
+    #     BW_WFS,
+    #     bw_json
+    # )
+    #
+    # with open(bw_json) as json_file:
+    #     data = json.load(json_file)
+    # xs, ys = [], []
+    # import matplotlib.pyplot as plt
+    # for feature in data["features"]:
+    #     x, y = feature['bbox'][:2]
+    #     date = feature['properties']['produktionsdatum'] or '1970-01-01'
+    #     xs.append(x)
+    #     ys.append(y)
+    #     # citygmls.append(
+    #     #     CityGMLWithDate(
+    #     #         url=feature["properties"]["xml"],
+    #     #         bundesland=NIEDERSACHSEN,
+    #     #         source_date=datetime.strptime(feature["properties"]["Aktualitaet"], "%Y-%m-%d %H:%M:%S"),
+    #     #         # coordinate_reference_system ?
+    #     #     )
+    #     # )
+    # # TODO: Show dates too?
+    # plt.scatter(xs, ys, s=0.1)
+    # plt.show()
+    from owslib.wfs import WebFeatureService
+    import geopandas as gpd
+    import io
+    from citygml_download import TMP_DIR
+
+# 1. Define the WFS URL and Feature Type
+    WFS_URL = "https://owsproxy.lgl-bw.de/owsproxy/wfs/WFS_LGL-BW_LoD2_Aktualitaet?"
+    FEATURE_TYPE = "verm:v_lod2_aktualitaet"
+
+# 2. Connect to the WFS
+    try:
+        # Use version 2.0.0 as it's the modern standard
+        wfs = WebFeatureService(url=WFS_URL, version="2.0.0")
+    except Exception as e:
+        print(f"Error connecting to WFS: {e}")
+        # Handle error or exit
+
+# 3. Request all features in GeoJSON format
+# We specify a large 'count' (MAXIMUM=1000000) to request all tiles.
+# Since the server handles streaming large files, this should work without manual pagination.
+    print(f"Requesting features for {FEATURE_TYPE}...")
+
+    try:
+        # Get the raw response content (GeoJSON string)
+        response = wfs.getfeature(
+            typename=FEATURE_TYPE,
+            outputFormat="application/json",
+            # count=1000000,  # Request a large number to ensure all tiles are returned
+        ).read()
+
+        # 4. Convert the GeoJSON response to a GeoPandas GeoDataFrame
+        # Use io.BytesIO to treat the response string as a file for geopandas
+        gdf = gpd.read_file(io.BytesIO(response), driver="GeoJSON")
+
+        # 5. Review and Save the Data
+
+        # Print the total number of tiles successfully loaded
+        print(f"\n✅ Download complete. Total tiles loaded: {len(gdf)}")
+
+        # Save the GeoDataFrame to a file
+        file_path = "bw_lod2_aktualitaet_tiles.geojson"
+        gdf.to_file(TMP_DIR / file_path, driver="GeoJSON")
+        print(f"Data saved to: {file_path}")
+
+        print(gdf.columns)
+
+        # Inspect the key columns
+        print("\nFirst 5 tiles and their Kachelnamen:")
+        print(gdf[["kachelname", "produktionsdatum", "dgm_datum"]].head())
+        print(gdf["dgm_datum"].describe())
+        print(gdf["produktionsdatum"].describe())
+        # Assuming your GeoDataFrame is named 'gdf'
+        import pandas as pd
+
+# 1. Ensure 'produktionsdatum' is in datetime format
+        gdf['produktionsdatum'] = pd.to_datetime(gdf['produktionsdatum'])
+
+# 2. Sort the GeoDataFrame by the more recent 'produktionsdatum'
+# This column effectively captures the latest update that led to file re-packaging.
+        most_recent_tiles = gdf.sort_values(by='produktionsdatum', ascending=False)
+
+# 3. Display the top 10 most recently updated tiles
+        print("Top 10 most recently updated tiles, sorted by 'produktionsdatum':")
+        print(most_recent_tiles[['kachelname', 'produktionsdatum', 'dgm_datum']].head(10))
+
+    except Exception as e:
+        print(f"\n❌ An error occurred during the GetFeature request or processing: {e}")
+        print("If the error is related to size, consider manually paginating the requests.")