adding sea polygons successful

plot.py

@@ -6,6 +6,34 @@ import geopandas as gpd
 import matplotlib.pyplot as plt
 from pathlib import Path
 
+def get_color(tags, geom_type):
+    """
+    Determine the color based on the tags and geometry type.
+    """
+    if isinstance(tags, dict):
+        if "natural" in tags:
+            if tags["natural"] == "water":
+                return "blue" if geom_type == "Polygon" else "darkblue"
+            if tags["natural"] in ["brush", "forest"]:
+                return "darkgreen"
+            if tags["natural"] == "coastline":
+                return "red"
+        if "landuse" in tags:
+            if tags["landuse"] == "farmland":
+                return "green"
+            if tags["landuse"] in ["residential", "commercial"]:
+                return "orange"
+            if tags["landuse"] == "industrial":
+                return "darkgrey"
+            return "yellow"
+        if "military" in tags:
+            return "lightred"
+        if "waterway" in tags:
+            return "darkblue"
+        if "railway" in tags:
+            return "black"
+    return "grey"  # Default color
+
 def plot_geometries(parquet_file: str, save_as: str):
     """
     Plots the geometries inside a GeoParquet file and saves the image instead of showing it.
@@ -25,9 +53,21 @@ def plot_geometries(parquet_file: str, save_as: str):
     if "geometry" not in gdf.columns:
         raise ValueError("❌ No 'geometry' column found in the file.")
     
-    # Plot the geometries
+    # Drop rows with empty geometries
+    gdf = gdf[gdf.geometry.notnull()]
+    
+    # Create a figure
     fig, ax = plt.subplots(figsize=(8, 8))
-    gdf.plot(ax=ax, edgecolor="black", facecolor="lightblue", alpha=0.5)
+    
+    # Classify colors efficiently
+    gdf["color"] = gdf.apply(lambda row: get_color(row.tags, row.geometry.geom_type), axis=1)
+    
+    # Separate by geometry types for efficient plotting
+    for geom_type, sub_gdf in gdf.groupby(gdf.geometry.geom_type):
+        sub_gdf.plot(ax=ax, edgecolor="black" if geom_type in ["Polygon", "MultiPolygon"] else None,
+                     facecolor=sub_gdf["color"] if geom_type in ["Polygon", "MultiPolygon"] else None,
+                     color=sub_gdf["color"] if geom_type in ["LineString", "MultiLineString"] else None,
+                     alpha=0.5, linewidth=2 if geom_type in ["LineString", "MultiLineString"] else None)
     
     # Add labels
     ax.set_title(f"Geometries in {parquet_file}")

water.py

@@ -3,143 +3,8 @@ import argparse
 import pandas as pd
 import geopandas as gpd
 import geohash
-import numpy as np
 import logging
-from shapely.geometry import box, LineString, MultiLineString, Point
+from shapely.geometry import box
-from shapely.ops import unary_union, linemerge, split
-import matplotlib.pyplot as plt
-
-def split_polygon_by_coastline(tile_polygon, coastline):
-    """
-    Splits the geohash tile polygon into two by the coastline and keeps the polygon on the 'right' side
-    based on the LineString's direction.
-    
-    Parameters:
-        - tile_polygon: The full geohash bounding box as a Polygon.
-        - coastline: The merged LineString coastline within the tile.
-    
-    Returns:
-        - The remaining water polygon after removing the 'left' side.
-    """
-    if coastline.is_empty or coastline.geom_type not in ["LineString", "MultiLineString"]:
-        return None  # No valid coastline
-
-    # Convert MultiLineString to the longest single LineString if necessary
-    if coastline.geom_type == "MultiLineString":
-        coastline = max(coastline.geoms, key=lambda g: g.length)
-
-    # Get the midpoint of the LineString
-    midpoint = coastline.interpolate(0.5, normalized=True)
-
-    # Compute the direction of the coastline (from first to last point)
-    coords = np.array(coastline.coords)
-    start, end = coords[0], coords[-1]
-    dx, dy = end[0] - start[0], end[1] - start[1]
-
-    # Compute normal vector (perpendicular to the coastline direction)
-    normal = np.array([-dy, dx])  # Rotate 90° counterclockwise
-    normal = normal / np.linalg.norm(normal)  # Normalize
-
-    # Compute a test point on the 'right' side
-    test_point = midpoint.x + normal[0], midpoint.y + normal[1]
-
-    # Perform the split
-    split_result = split(tile_polygon, coastline)
-
-    if not split_result or len(split_result.geoms) < 2:
-        return None  # No valid split
-
-    # Choose the polygon that contains the test point
-    for poly in split_result.geoms:
-        if poly.contains(Point(test_point)):
-            return poly  # Keep the 'right' side
-
-    return None  # Fail-safe
-
-
-def group_continuous_coastline(segments, gap_threshold=0.001):
-    """
-    Groups and merges only those LineStrings that are properly connected to each other.
-    Separate coastlines remain unmerged (e.g., islands).
-    
-    Parameters:
-        - segments: List of LineStrings.
-        - gap_threshold: Max distance (degrees) between segment endpoints to be considered connected.
-    
-    Returns:
-        - List of merged continuous LineStrings.
-    """
-    connected_groups = []
-    remaining_segments = list(segments)
-
-    while remaining_segments:
-        # Start a new group with one segment
-        group = [remaining_segments.pop(0)]
-
-        added = True
-        while added:
-            added = False
-            for seg in remaining_segments[:]:
-                if is_connected_to_group(seg, group, gap_threshold):
-                    group.append(seg)
-                    remaining_segments.remove(seg)
-                    added = True  # Continue expanding the group
-
-        # Merge each connected group separately
-        if len(group) > 1:
-            merged = linemerge(MultiLineString(group))
-        else:
-            merged = group[0]
-
-        connected_groups.append(merged)
-
-    return connected_groups
-
-
-def is_connected_to_group(line, group, gap_threshold=0.001):
-    """
-    Checks if a LineString is connected to any other in a group.
-    A connection means they share an endpoint within a small gap threshold.
-    """
-    line_start, line_end = Point(line.coords[0]), Point(line.coords[-1])
-
-    for other in group:
-        other_start, other_end = Point(other.coords[0]), Point(other.coords[-1])
-        if (
-            line_start.distance(other_end) < gap_threshold or
-            line_end.distance(other_start) < gap_threshold
-        ):
-            return True  # They are connected within the gap threshold
-
-    return False
-
-def plot_coastline_segments(geohash_code, fragmented_coastline, output_folder="debug_plots"):
-    """Saves a PNG plot of coastline segments for debugging."""
-    os.makedirs(output_folder, exist_ok=True)  # Ensure the output folder exists
-    output_file = os.path.join(output_folder, f"{geohash_code}.png")
-
-    fig, ax = plt.subplots(figsize=(6, 6))
-    colors = ['red', 'blue', 'green', 'purple', 'orange', 'cyan']
-
-    # Convert to list if MultiLineString
-    if fragmented_coastline.geom_type == "MultiLineString":
-        segments = list(fragmented_coastline.geoms)
-    else:
-        segments = [fragmented_coastline]
-
-    for idx, segment in enumerate(segments):
-        x, y = segment.xy
-        ax.plot(x, y, color=colors[idx % len(colors)], linewidth=2, label=f"Segment {idx+1}")
-
-    ax.set_title(f"Coastline Segments in {geohash_code}")
-    ax.set_xlabel("Longitude")
-    ax.set_ylabel("Latitude")
-    ax.legend()
-    
-    plt.savefig(output_file, dpi=300, bbox_inches="tight")  # Save as PNG
-    plt.close(fig)  # Close the plot to free memory
-
-    print(f"✅ Saved coastline debug plot: {output_file}")
 
 # Configure logging
 logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.INFO)
@@ -163,52 +28,17 @@ def geohash_bbox(geohash_code):
 def ensure_crs_consistency(gdf):
     """Ensures the GeoDataFrame is in EPSG:4326 to prevent CRS mismatches."""
     if gdf.crs is None or gdf.crs.to_string() != "EPSG:4326":
-        logger.info(f"Converting CRS of {gdf} to EPSG:4326")
         gdf = gdf.to_crs("EPSG:4326")
     return gdf
 
-def check_coastline_continuity(coastline_gdf, geohash_code):
+def process_geohash_files(water_gdf, geohash_files):
-    """Checks if the coastline forms a continuous LineString and attempts to merge fragmented parts."""
+    """Processes each geohash file: adds water polygons from water.parquet."""
-    if coastline_gdf.empty:
+    # Precompute bounding boxes for all geohash codes
-        logger.warning(f"Geohash {geohash_code}: No coastline to check.")
+    geohash_bboxes = {geohash_code: geohash_bbox(geohash_code) for geohash_code in geohash_files.keys()}
-        return False
 
-    # Convert to a list of valid LineStrings
+    # Spatial index for faster spatial queries
-    coastline_list = [geom for geom in coastline_gdf.geometry if geom and geom.geom_type in ["LineString", "MultiLineString"]]
+    water_sindex = water_gdf.sindex
 
-    if not coastline_list:
-        logger.error(f"Geohash {geohash_code}: No valid LineString geometries found.")
-        return False
-
-    # Merge the coastline segments (ensuring not empty)
-    merged_coastlines = group_continuous_coastline(coastline_list)
-
-    if not merged_coastlines:
-        logger.error(f"Geohash {geohash_code}: Merging failed. No valid coastlines found.")
-        merged_coastline = MultiLineString([])  # Return an empty MultiLineString
-    else:
-        merged_coastline = max(merged_coastlines, key=lambda ls: ls.length)
-
-    # Ensure the result is valid
-    if merged_coastline.is_empty:
-        logger.error(f"Geohash {geohash_code}: Merged coastline is empty.")
-        return False
-
-    if isinstance(merged_coastline, LineString):
-        logger.info(f"Geohash {geohash_code}: Coastline is continuous.")
-        plot_coastline_segments(geohash_code, merged_coastline)
-        return True
-
-    if isinstance(merged_coastline, MultiLineString):
-        logger.warning(f"Geohash {geohash_code}: Coastline remains fragmented ({len(merged_coastline.geoms)} parts).")
-        plot_coastline_segments(geohash_code, merged_coastline)
-        return False
-
-    logger.error(f"Geohash {geohash_code}: Unexpected coastline geometry ({merged_coastline.geom_type}).")
-    return False
-
-def process_geohash_files(coastline_gdf, geohash_files):
-    """Processes each geohash file: updates coastline and adds water polygons."""
     for geohash_code, file_path in geohash_files.items():
         logger.info(f"Processing {file_path}...")
 
@@ -221,65 +51,48 @@ def process_geohash_files(coastline_gdf, geohash_files):
 
         existing_gdf = ensure_crs_consistency(existing_gdf)
 
-        # Remove existing coastline features
+        # Get the precomputed bounding box for this geohash
-        filtered_gdf = existing_gdf[~((existing_gdf.geometry.type.isin(["LineString", "MultiLineString"])) &
+        bbox = geohash_bboxes[geohash_code]
-                                      (existing_gdf["tags"].apply(lambda tags: isinstance(tags, dict) and tags.get("natural") == "coastline")))]
 
-        # Get coastline geometries for this geohash
+        # Use spatial index to find intersecting water polygons
-        bbox = geohash_bbox(geohash_code)
+        possible_matches_index = list(water_sindex.intersection(bbox.bounds))
-        new_coastline = coastline_gdf[coastline_gdf.intersects(bbox)].copy()
+        possible_matches = water_gdf.iloc[possible_matches_index]
+        water_in_tile = possible_matches[possible_matches.intersects(bbox)].copy()
 
-        if new_coastline.empty:
+        if water_in_tile.empty:
-            logger.info(f"No coastline found for {geohash_code}. Skipping coastline update.")
+            logger.info(f"No water found for {geohash_code}. Skipping water addition.")
             continue
 
-        # Merge continuous coastline segments
+        # Clip water polygons to the geohash tile boundary
-        merged_coastlines = group_continuous_coastline(new_coastline.geometry)
+        water_in_tile["geometry"] = water_in_tile.intersection(bbox)
-        if merged_coastlines:
-            merged_coastline = max(merged_coastlines, key=lambda ls: ls.length)  # Pick the longest coastline
-        else:
-            merged_coastline = None
 
-        if merged_coastline is None:
+        # Assign "water" metadata
-            logger.warning(f"Geohash {geohash_code}: Coastline merging failed.")
+        water_in_tile["tags"] = [{"natural": "water", "water": "sea"}] * len(water_in_tile)
-            continue
 
-        logger.info(f"Adding water polygon for {geohash_code}...")
+        # Ensure columns match before merging
+        for col in ["feature_id", "tags"]:
+            if col not in existing_gdf.columns:
+                existing_gdf[col] = None
+            if col not in water_in_tile.columns:
+                water_in_tile[col] = None
 
-        # Create a full water polygon covering the entire geohash tile
+        # Merge updated water polygons into geohash file
-        tile_polygon = bbox
+        updated_gdf = gpd.GeoDataFrame(pd.concat([existing_gdf, water_in_tile], ignore_index=True), crs="EPSG:4326")
-
-        # Split the tile polygon using the coastline
-        water_polygon = split_polygon_by_coastline(tile_polygon, merged_coastline)
-
-        if water_polygon:
-            # Assign water tags
-            water_gdf = gpd.GeoDataFrame(
-                {"geometry": [water_polygon], "tags": [{"natural": "water", "water": "sea"}]},
-                crs="EPSG:4326"
-            )
-
-            # Append the water polygon to the updated geohash file
-            updated_gdf = gpd.GeoDataFrame(pd.concat([filtered_gdf, new_coastline, water_gdf], ignore_index=True), crs="EPSG:4326")
-            logger.info(f"Water polygon added for {geohash_code}.")
-        else:
-            logger.warning(f"Could not split water polygon for {geohash_code}. Skipping water addition.")
-            updated_gdf = gpd.GeoDataFrame(pd.concat([filtered_gdf, new_coastline], ignore_index=True), crs="EPSG:4326")
 
         # Save back to parquet
         updated_gdf.to_parquet(file_path, index=False)
-        logger.info(f"Updated geohash file {file_path}.")
+        logger.info(f"Updated geohash file {file_path} with water polygons.")
 
 def main():
-    parser = argparse.ArgumentParser(description="Update coastline geometries in existing geohash parquet files.")
+    parser = argparse.ArgumentParser(description="Add water polygons to geohash parquet files from a global water dataset.")
-    parser.add_argument("coastline", help="Path to the coastline.parquet file.")
+    parser.add_argument("water", help="Path to the water.parquet file.")
     parser.add_argument("geohash_folder", help="Path to the folder containing <geohash>.parquet files.")
     
     args = parser.parse_args()
 
-    # Load coastline dataset
+    # Load water dataset
-    logger.info(f"Loading coastline dataset from {args.coastline}...")
+    logger.info(f"Loading water dataset from {args.water}...")
-    coastline_gdf = ensure_crs_consistency(gpd.read_parquet(args.coastline))
+    water_gdf = ensure_crs_consistency(gpd.read_parquet(args.water))
 
     # Get existing geohash parquet files
     geohash_files = load_geohash_files(args.geohash_folder)
@@ -289,9 +102,9 @@ def main():
         return
 
     # Process geohash files
-    process_geohash_files(coastline_gdf, geohash_files)
+    process_geohash_files(water_gdf, geohash_files)
 
-    logger.info("Processing complete. All coastline updates applied.")
+    logger.info("Processing complete. All water polygons added.")
 
 if __name__ == "__main__":
     main()