From f2931af8dd11f2b469c216459caf911831b2c492 Mon Sep 17 00:00:00 2001
From: randogoth <kflux@posteo.de>
Date: Wed, 26 Feb 2025 14:01:09 +0000
Subject: [PATCH] slice retains tags, lookup function

---
 import.py   | 10 ++++---
 install.txt |  2 +-
 lookup.py   | 81 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 slice.py    |  2 +-
 4 files changed, 89 insertions(+), 6 deletions(-)
 create mode 100644 lookup.py

diff --git a/import.py b/import.py
index ff8fcf3..d2b6c81 100644
--- a/import.py
+++ b/import.py
@@ -86,7 +86,9 @@ def bbox_geohashes(bbox: Tuple[float, float, float, float]) -> List[str]:
 
     return sorted(geohashes)
 
-remove_points('compact.parquet', 'processed.parquet')
-bbox = bbox('israel.pbf')
-hashes = bbox_geohashes(bbox)
-print(hashes)
\ No newline at end of file
+# remove_points('compact.parquet', 'processed.parquet')
+# bbox = bbox('israel.pbf')
+# hashes = bbox_geohashes(bbox)
+# print(hashes)
+
+check_data("geohash/sv2c.parquet")
\ No newline at end of file
diff --git a/install.txt b/install.txt
index 30375af..976371e 100644
--- a/install.txt
+++ b/install.txt
@@ -1,4 +1,4 @@
 ## packages
 sudo dnf install osmium-tool jq python3-devel pip
 
-pip install quackosm[cli] psycopg2-binary pyosmium geoalchemy2 geojson shapely
\ No newline at end of file
+pip install quackosm[cli] psycopg2-binary pyosmium geoalchemy2 geojson shapely pygeohash pyarrow pandas click
diff --git a/lookup.py b/lookup.py
new file mode 100644
index 0000000..1e3d408
--- /dev/null
+++ b/lookup.py
@@ -0,0 +1,81 @@
+import json
+import click
+import pygeohash as pgh
+import pyarrow.parquet as pq
+from shapely.geometry import Point, LineString, MultiLineString
+from shapely import wkb
+from pyproj import Geod
+
+@click.command()
+@click.argument("latitude", type=float)
+@click.argument("longitude", type=float)
+@click.option("--distance", type=float, default=10, help="Distance threshold in meters for nearby lines.")
+@click.option("--parquet-dir", type=str, default="geohash", help="Directory containing geohash Parquet files.")
+def geohash_lookup(latitude, longitude, distance, parquet_dir):
+    """
+    Finds polygons containing the given point and lines within a threshold (in meters).
+    """
+
+    # Compute geohash (precision can be adjusted)
+    geohash = pgh.encode(latitude, longitude, precision=4)  # Adjust precision as needed
+    parquet_file = f"{parquet_dir}/{geohash}.parquet"
+
+    try:
+        # Load the Parquet file
+        table = pq.read_table(parquet_file)
+        df = table.to_pandas()
+    except FileNotFoundError:
+        click.echo(json.dumps({"error": f"No Parquet file found for geohash {geohash}."}))
+        return
+
+    point = Point(longitude, latitude)
+    geod = Geod(ellps="WGS84")  # Accurate geodetic distance calculations
+
+    found_geometries = []
+
+    for _, row in df.iterrows():
+        # Decode geometry from WKB
+        try:
+            geom = wkb.loads(bytes(row["geometry"]))  # Ensure WKB is properly read
+        except Exception as e:
+            click.echo(json.dumps({"error": f"Failed to decode geometry: {str(e)}"}))
+            continue
+
+        # ✅ Ensure `tags` is always a dictionary
+        metadata = {}
+        if isinstance(row["tags"], dict):
+            metadata = {k: v for k, v in row["tags"].items() if v is not None}  # Remove empty values
+
+        # ✅ Check if the point is inside a polygon
+        if geom.contains(point):
+            found_geometries.append({"type": "polygon", **metadata})
+
+        # ✅ Handle LineString and MultiLineString
+        elif isinstance(geom, LineString):
+            try:
+                min_distance = min(geod.inv(point.x, point.y, p[0], p[1])[2] for p in geom.coords)
+                if min_distance <= distance:
+                    found_geometries.append({"type": "line", **metadata})
+            except Exception as e:
+                click.echo(json.dumps({"error": f"Failed distance calculation (LineString): {str(e)}"}))
+                continue
+
+        elif isinstance(geom, MultiLineString):
+            try:
+                # ✅ Iterate through all LineStrings inside the MultiLineString
+                min_distance = min(
+                    min(geod.inv(point.x, point.y, p[0], p[1])[2] for p in line.coords)
+                    for line in geom.geoms
+                )
+
+                if min_distance <= distance:
+                    found_geometries.append({"type": "multi_line", **metadata})
+            except Exception as e:
+                click.echo(json.dumps({"error": f"Failed distance calculation (MultiLineString): {str(e)}"}))
+                continue
+
+    # Output JSON result
+    click.echo(json.dumps(found_geometries, indent=2, ensure_ascii=False))  # Ensure proper Unicode display
+
+if __name__ == "__main__":
+    geohash_lookup()
diff --git a/slice.py b/slice.py
index 0a9e66d..dce9115 100644
--- a/slice.py
+++ b/slice.py
@@ -97,7 +97,7 @@ def slice_and_split_geoparquet(input_file: str, output_dir: str):
 
         # Clip geometries inside this geohash box
         filtered_data = con.execute(f"""
-            SELECT feature_id, rowid, 
+            SELECT feature_id, tags, 
                    ST_AsWKB(ST_Intersection(geometry, 
                                             ST_MakeEnvelope({geohash_bbox["w"]}, {geohash_bbox["s"]}, 
                                                             {geohash_bbox["e"]}, {geohash_bbox["n"]})))