slice retains tags, lookup function

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')
+# remove_points('compact.parquet', 'processed.parquet')
-bbox = bbox('israel.pbf')
+# bbox = bbox('israel.pbf')
-hashes = bbox_geohashes(bbox)
+# hashes = bbox_geohashes(bbox)
-print(hashes)
+# print(hashes)
+
+check_data("geohash/sv2c.parquet")

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
+pip install quackosm[cli] psycopg2-binary pyosmium geoalchemy2 geojson shapely pygeohash pyarrow pandas click

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()

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"]})))