init

.gitignore

@@ -0,0 +1,2 @@
+*.parquet
+*.pbf

import.py

@@ -0,0 +1,92 @@
+import duckdb
+from pathlib import Path
+import geohash
+import osmium
+import subprocess
+from pprint import pprint
+from typing import List, Tuple, Dict
+
+def remove_points(input_file: str, output_file: str):
+    try:
+        input_path, output_path = Path(input_file), Path(output_file)
+        if not input_path.exists():
+            raise FileNotFoundError(f"File '{input_file}' not found.")
+
+        con = duckdb.connect()
+        con.execute("INSTALL spatial; LOAD spatial;")
+
+        con.execute(f"""
+            CREATE VIEW geoparquet_filtered AS 
+            SELECT * FROM read_parquet('{input_path}')
+            WHERE ST_GeometryType(geometry) IS NOT NULL
+              AND ST_GeometryType(geometry) != 'POINT'
+        """)
+
+        if con.execute("SELECT COUNT(*) FROM geoparquet_filtered").fetchone()[0] == 0:
+            raise ValueError("No non-POINT geometries found.")
+
+        output_path.parent.mkdir(parents=True, exist_ok=True)
+        con.execute(f"COPY geoparquet_filtered TO '{output_path}' (FORMAT PARQUET)")
+        print(f"✅ Saved to {output_path}")
+
+    except Exception as e:
+        print(f"❌ Error: {e}")
+
+    finally:
+        if 'con' in locals():
+            con.close()
+
+def bbox(file_path: str) -> Dict[str, float]:
+    try:
+        result = subprocess.run(
+            ["osmium", "fileinfo", "-e", file_path],
+            capture_output=True, text=True, check=True
+        )
+        for line in result.stdout.split("\n"):
+            if "Bounding box:" in line:
+                bbox_str = line.strip().replace("Bounding box:", "").strip()
+                # Convert to a tuple of floats
+                bbox_tuple = tuple(map(float, bbox_str.strip("()").split(",")))
+                return bbox_tuple
+    except subprocess.CalledProcessError as e:
+        print(f"❌ Error running osmium: {e}")
+
+def check_data(file: str):
+    con = duckdb.connect()
+    con.execute("INSTALL spatial; LOAD spatial;")
+    con.execute(f"CREATE VIEW geoparquet AS SELECT * FROM read_parquet('{file}')")
+    columns = con.execute("DESCRIBE geoparquet").fetchdf()
+    print("Columns in the file:\n", columns)
+    dc = con.execute("SELECT count(*) FROM geoparquet LIMIT 10").fetchdf()
+    print("\nTotal rows:", dc)
+    df = con.execute("SELECT * FROM geoparquet LIMIT 10").fetchdf()
+    print("\nFirst 10 rows:\n", df)
+    try:
+        geo_types = con.execute("SELECT DISTINCT ST_GeometryType(geometry) AS geom_type FROM geoparquet").fetchdf()
+        print("\nDetected Geometry Types:\n", geo_types)
+    except Exception as e:
+        print("\nCould not determine geometry type. Possible reason: No valid geometry column found.\n", e)
+    con.close()
+
+def bbox_geohashes(bbox: Tuple[float, float, float, float]) -> List[str]:
+    min_x, min_y, max_x, max_y = bbox
+    geohashes = set()
+
+    # Define step size based on 4-character Geohash precision (~39km × 19.5km)
+    step_size = 0.25  # Approximate step size in degrees (varies slightly by latitude)
+
+    lat = min_y
+    while lat <= max_y:
+        lon = min_x
+        while lon <= max_x:
+            gh = geohash.encode(lat, lon, precision=4)
+            geohashes.add(gh)
+            lon += step_size  # Move eastward
+        lat += step_size  # Move northward
+
+    return sorted(geohashes)
+
+remove_points('compact.parquet', 'processed.parquet')
+bbox = bbox('israel.pbf')
+hashes = bbox_geohashes(bbox)
+print(hashes)

import_osm.py

@@ -0,0 +1,58 @@
+import json
+import psycopg2
+from shapely.geometry import shape
+
+def import_geojson(filepath):
+    conn = psycopg2.connect(
+        dbname="gis_data",
+        user="osm",
+        password="osm",
+        host="localhost"
+    )
+    cur = conn.cursor()
+
+    try:
+        with open(filepath) as f:
+            data = json.load(f)
+
+        if data["type"] != "FeatureCollection":
+            raise ValueError("GeoJSON must be a FeatureCollection")
+
+        for idx, feature in enumerate(data["features"]):
+            try:
+                geom_type = feature["geometry"]["type"]
+                properties = feature["properties"]
+                geom = shape(feature["geometry"])
+                wkt = geom.wkt
+
+                # Insert lines
+                if geom_type == "LineString":
+                    cur.execute("""
+                        INSERT INTO osm_lines (geometry, tags)
+                        VALUES (ST_GeomFromText(%s, 4326), %s)
+                    """, (wkt, json.dumps(properties)))
+
+                # Insert polygons AND multipolygons
+                elif geom_type in ("Polygon", "MultiPolygon"):
+                    cur.execute("""
+                        INSERT INTO osm_polygons (geometry, tags)
+                        VALUES (ST_GeomFromText(%s, 4326), %s)
+                    """, (wkt, json.dumps(properties)))
+
+                else:
+                    print(f"Skipping unsupported geometry type: {geom_type}")
+
+            except Exception as e:
+                print(f"Error processing feature {idx}: {e}")
+                continue
+
+        conn.commit()
+
+    except json.JSONDecodeError:
+        print("Invalid JSON file. Ensure filtered.geojson is correctly generated.")
+    finally:
+        cur.close()
+        conn.close()
+
+if __name__ == "__main__":
+    import_geojson("israel.geojson")

install.txt

@@ -0,0 +1,30 @@
+## packages
+sudo dnf install postgresql-server postgresql-contrib postgis osmium-tool jq python3-devel pip
+
+pip install quackosm[cli] psycopg2-binary pyosmium geoalchemy2 geojson shapely
+
+## PostgreSQL setup
+sudo postgresql-setup --initdb
+sudo -u postgres psql
+
+CREATE USER osm WITH PASSWORD 'osm';
+CREATE DATABASE gis_data OWNER osm;
+\q
+sudo -u postgres psql -d gis_data
+
+CREATE TABLE osm_lines (
+    id SERIAL PRIMARY KEY,
+    geometry GEOMETRY(LINESTRING, 4326),
+    tags JSONB
+);
+
+CREATE TABLE osm_polygons (
+    id SERIAL PRIMARY KEY,
+    geometry GEOMETRY(POLYGON, 4326),
+    tags JSONB
+);
+
+CREATE INDEX idx_osm_lines_geometry ON osm_lines USING GIST(geometry);
+CREATE INDEX idx_osm_polygons_geometry ON osm_polygons USING GIST(geometry);
+CREATE INDEX idx_osm_lines_tags ON osm_lines USING GIN(tags);
+CREATE INDEX idx_osm_polygons_tags ON osm_polygons USING GIN(tags);

plot.py

@@ -0,0 +1,54 @@
+import argparse
+import matplotlib
+matplotlib.use("Agg")  # Use Agg backend for non-GUI environments
+
+import geopandas as gpd
+import matplotlib.pyplot as plt
+from pathlib import Path
+
+def plot_geometries(parquet_file: str, save_as: str):
+    """
+    Plots the geometries inside a GeoParquet file and saves the image instead of showing it.
+
+    Args:
+        parquet_file (str): Path to the Parquet file containing geometries.
+        save_as (str): Path to save the output plot image.
+    """
+    parquet_path = Path(parquet_file)
+    if not parquet_path.exists():
+        raise FileNotFoundError(f"❌ File '{parquet_file}' not found.")
+
+    # Load the Parquet file into a GeoDataFrame
+    gdf = gpd.read_parquet(parquet_path)
+
+    # Ensure it has a valid geometry column
+    if "geometry" not in gdf.columns:
+        raise ValueError("❌ No 'geometry' column found in the file.")
+
+    # Plot the geometries
+    fig, ax = plt.subplots(figsize=(8, 8))
+    gdf.plot(ax=ax, edgecolor="black", facecolor="lightblue", alpha=0.5)
+
+    # Add labels
+    ax.set_title(f"Geometries in {parquet_file}")
+    ax.set_xlabel("Longitude")
+    ax.set_ylabel("Latitude")
+
+    # Save the figure instead of showing it
+    plt.savefig(save_as, dpi=300)
+    print(f"✅ Plot saved as {save_as}")
+
+def main():
+    parser = argparse.ArgumentParser(description="Plot geometries from a GeoParquet file and save as an image.")
+    parser.add_argument("-i", "--input", required=True, help="Path to the input GeoParquet file.")
+    parser.add_argument("-o", "--output", required=True, help="Path to save the output plot image.")
+    
+    args = parser.parse_args()
+
+    try:
+        plot_geometries(args.input, args.output)
+    except Exception as e:
+        print(f"❌ Error: {e}")
+
+if __name__ == "__main__":
+    main()

slice.py

@@ -0,0 +1,155 @@
+import argparse
+import duckdb
+import geohash
+import geopandas as gpd
+import pandas as pd
+import pyarrow.parquet as pq
+import pyarrow as pa
+from pathlib import Path
+from shapely import wkb
+from typing import List, Dict
+
+def get_geohashes_from_bbox(min_x, min_y, max_x, max_y) -> List[str]:
+    """Generates all 4-character Geohashes that intersect a bounding box."""
+    geohashes = set()
+    step_size = 0.25  # 4-character Geohash grid size
+
+    lat = min_y
+    while lat <= max_y:
+        lon = min_x
+        while lon <= max_x:
+            gh = geohash.encode(lat, lon, precision=4)
+            geohashes.add(gh)
+            lon += step_size
+        lat += step_size
+
+    return sorted(geohashes)
+
+def append_to_parquet(file_path: Path, new_data: gpd.GeoDataFrame):
+    """Appends new data to an existing Parquet file while ensuring correct merging."""
+    if file_path.exists():
+        try:
+            # ✅ Load existing data
+            existing_data = gpd.read_parquet(file_path)
+            
+            # ✅ Merge with new data
+            combined_data = pd.concat([existing_data, new_data], ignore_index=True)
+
+            # ✅ Drop duplicates *only if feature_id exists*
+            if "feature_id" in combined_data.columns:
+                combined_data = combined_data.drop_duplicates(subset="feature_id", keep="last")
+
+            # ✅ Reset index before saving
+            combined_data = combined_data.reset_index(drop=True)
+
+            # ✅ Save back to the same file
+            combined_data.to_parquet(file_path, index=False)
+            print(f"✅ Appended new data to {file_path}")
+
+        except Exception as e:
+            print(f"❌ Error while merging {file_path}: {e}")
+
+    else:
+        # If the file doesn't exist, create it
+        new_data.to_parquet(file_path, index=False)
+        print(f"✅ Created new file: {file_path}")
+
+def slice_and_split_geoparquet(input_file: str, output_dir: str):
+    """
+    Splits the GeoParquet file into multiple files based on 4-character Geohash boxes.
+    Appends new data to existing Geohash Parquet files while preventing duplicate geometries.
+    """
+    input_path, output_path = Path(input_file), Path(output_dir)
+    if not input_path.exists():
+        raise FileNotFoundError(f"File '{input_file}' not found.")
+    output_path.mkdir(parents=True, exist_ok=True)
+
+    print(f"📂 Processing '{input_file}', output will be saved in '{output_dir}'...")
+
+    con = duckdb.connect()
+    con.execute("INSTALL spatial; LOAD spatial;")
+
+    # Create a temporary table of all non POINT geometries
+    con.execute(f"CREATE TEMP TABLE geoparquet AS SELECT * FROM read_parquet('{input_path}') WHERE ST_GeometryType(geometry) IS NOT NULL AND ST_GeometryType(geometry) != 'POINT'")
+
+    # Compute bounding box for each geometry
+    bbox_results = con.execute("""
+        SELECT feature_id, rowid, ST_XMin(ST_Envelope(geometry)), ST_YMin(ST_Envelope(geometry)),
+                     ST_XMax(ST_Envelope(geometry)), ST_YMax(ST_Envelope(geometry))
+        FROM geoparquet;
+    """).fetchall()
+
+    # Compute all intersecting geohashes
+    geohash_mapping: Dict[str, List[int]] = {}
+
+    for feature_id, rowid, min_x, min_y, max_x, max_y in bbox_results:
+        intersecting_geohashes = get_geohashes_from_bbox(min_x, min_y, max_x, max_y)
+        for geohash_code in intersecting_geohashes:
+            if geohash_code not in geohash_mapping:
+                geohash_mapping[geohash_code] = []
+            geohash_mapping[geohash_code].append(rowid)
+
+    # Dictionary to store in-memory results before writing
+    geohash_data: Dict[str, gpd.GeoDataFrame] = {}
+
+    for geohash_code, rowids in geohash_mapping.items():
+        geohash_bbox = geohash.bbox(geohash_code)
+
+        # Clip geometries inside this geohash box
+        filtered_data = con.execute(f"""
+            SELECT feature_id, rowid, 
+                   ST_AsWKB(ST_Intersection(geometry, 
+                                            ST_MakeEnvelope({geohash_bbox["w"]}, {geohash_bbox["s"]}, 
+                                                            {geohash_bbox["e"]}, {geohash_bbox["n"]}))) 
+                   AS clipped_geom 
+            FROM geoparquet 
+            WHERE rowid IN ({','.join(map(str, rowids))})
+            AND ST_Intersects(geometry, 
+                              ST_MakeEnvelope({geohash_bbox["w"]}, {geohash_bbox["s"]}, 
+                                              {geohash_bbox["e"]}, {geohash_bbox["n"]}));
+        """).fetchdf()
+
+        if not filtered_data.empty:
+            # ✅ Remove non-WKB values
+            filtered_data = filtered_data[filtered_data["clipped_geom"].apply(lambda x: isinstance(x, (bytes, bytearray)))]
+
+            if not filtered_data.empty:
+                # ✅ Convert WKB to Shapely geometries safely
+                filtered_data["geometry"] = filtered_data["clipped_geom"].apply(
+                    lambda x: wkb.loads(bytes(x)) if isinstance(x, (bytes, bytearray)) else None
+                )
+                filtered_data.drop(columns=["clipped_geom"], inplace=True)
+
+                # Remove rows where geometry conversion failed
+                filtered_data = filtered_data.dropna(subset=["geometry"])
+
+                if not filtered_data.empty:
+                    # Convert to GeoPandas GeoDataFrame
+                    gdf = gpd.GeoDataFrame(filtered_data, geometry="geometry", crs="EPSG:4326")
+
+                    # Store in memory
+                    geohash_data[geohash_code] = gdf
+
+    # Step 4: Write each geohash's data **only once** and prevent duplicates
+    for geohash_code, gdf in geohash_data.items():
+        geohash_file = output_path / f"{geohash_code}.parquet"
+        append_to_parquet(geohash_file, gdf)
+        print(f"✅ Updated: {geohash_file}")
+
+    con.close()
+    print("🎉 Processing complete!")
+
+def main():
+    parser = argparse.ArgumentParser(description="Split a GeoParquet file into Geohash tiles.")
+    parser.add_argument("-i", "--input", required=True, help="Path to the input GeoParquet file.")
+    parser.add_argument("-o", "--output", required=True, help="Directory to save output Geohash tiles.")
+    
+    args = parser.parse_args()
+
+    try:
+        slice_and_split_geoparquet(args.input, args.output)
+    except Exception as e:
+        print(f"❌ Error: {e}")
+
+if __name__ == "__main__":
+    main()

tags.json

@@ -0,0 +1,19 @@
+{
+    "name:en": true,
+    "name": true,
+    "aeroway": true,
+    "boundary": ["aboriginal_lands", "border_zone", "forest", "hazard", "national_park", "protected_area", "disputed"],
+    "geological": true,
+    "highway": ["motorway", "trunk", "primary", "raceway"],
+    "historic": true,
+    "landuse": true,
+    "leisure": ["dog_park", "garden", "miniature_golf", "nature_reserve", "park", "playground", "stadium", "swimming_pool", "water_park"],
+    "military": ["airfield", "base", "bunker", "barracks", "danger_area", "nuclear_explosion_site", "range", "training_area"],
+    "natural": true,
+    "place": ["plot", "farm", "island", "islet", "sea", "ocean"],
+    "power": ["plant"],
+    "railway": ["abandoned", "construction", "disused", "rail"],
+    "tourism": ["attraction", "yes"],
+    "water": true,
+    "waterway": true
+}