From cdabeb730bd4e5d032184a39e2f6ee972964a18e Mon Sep 17 00:00:00 2001 From: randogoth Date: Tue, 4 Mar 2025 15:06:11 +0000 Subject: [PATCH] adding sea polygons successful --- plot.py | 52 +++++++++-- water.py | 259 ++++++++----------------------------------------------- 2 files changed, 82 insertions(+), 229 deletions(-) diff --git a/plot.py b/plot.py index f306067..8e44434 100644 --- a/plot.py +++ b/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. @@ -24,16 +52,28 @@ def plot_geometries(parquet_file: str, save_as: str): # 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 + + # 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}") 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}") @@ -51,4 +91,4 @@ def main(): print(f"❌ Error: {e}") if __name__ == "__main__": - main() + main() \ No newline at end of file diff --git a/water.py b/water.py index 027eb5e..49533d3 100644 --- a/water.py +++ b/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.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}") +from shapely.geometry import box # 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): - """Checks if the coastline forms a continuous LineString and attempts to merge fragmented parts.""" - if coastline_gdf.empty: - logger.warning(f"Geohash {geohash_code}: No coastline to check.") - return False +def process_geohash_files(water_gdf, geohash_files): + """Processes each geohash file: adds water polygons from water.parquet.""" + # Precompute bounding boxes for all geohash codes + geohash_bboxes = {geohash_code: geohash_bbox(geohash_code) for geohash_code in geohash_files.keys()} - # Convert to a list of valid LineStrings - coastline_list = [geom for geom in coastline_gdf.geometry if geom and geom.geom_type in ["LineString", "MultiLineString"]] + # Spatial index for faster spatial queries + 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 - filtered_gdf = existing_gdf[~((existing_gdf.geometry.type.isin(["LineString", "MultiLineString"])) & - (existing_gdf["tags"].apply(lambda tags: isinstance(tags, dict) and tags.get("natural") == "coastline")))] + # Get the precomputed bounding box for this geohash + bbox = geohash_bboxes[geohash_code] - # Get coastline geometries for this geohash - bbox = geohash_bbox(geohash_code) - new_coastline = coastline_gdf[coastline_gdf.intersects(bbox)].copy() + # Use spatial index to find intersecting water polygons + possible_matches_index = list(water_sindex.intersection(bbox.bounds)) + possible_matches = water_gdf.iloc[possible_matches_index] + water_in_tile = possible_matches[possible_matches.intersects(bbox)].copy() - if new_coastline.empty: - logger.info(f"No coastline found for {geohash_code}. Skipping coastline update.") + if water_in_tile.empty: + logger.info(f"No water found for {geohash_code}. Skipping water addition.") continue - # Merge continuous coastline segments - merged_coastlines = group_continuous_coastline(new_coastline.geometry) - if merged_coastlines: - merged_coastline = max(merged_coastlines, key=lambda ls: ls.length) # Pick the longest coastline - else: - merged_coastline = None + # Clip water polygons to the geohash tile boundary + water_in_tile["geometry"] = water_in_tile.intersection(bbox) - if merged_coastline is None: - logger.warning(f"Geohash {geohash_code}: Coastline merging failed.") - continue + # Assign "water" metadata + water_in_tile["tags"] = [{"natural": "water", "water": "sea"}] * len(water_in_tile) - 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 - tile_polygon = bbox - - # 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") + # Merge updated water polygons into geohash file + updated_gdf = gpd.GeoDataFrame(pd.concat([existing_gdf, water_in_tile], 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.add_argument("coastline", help="Path to the coastline.parquet file.") + parser = argparse.ArgumentParser(description="Add water polygons to geohash parquet files from a global water dataset.") + parser.add_argument("water", help="Path to the water.parquet file.") parser.add_argument("geohash_folder", help="Path to the folder containing .parquet files.") args = parser.parse_args() - # Load coastline dataset - logger.info(f"Loading coastline dataset from {args.coastline}...") - coastline_gdf = ensure_crs_consistency(gpd.read_parquet(args.coastline)) + # Load water dataset + logger.info(f"Loading water dataset from {args.water}...") + 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() + main() \ No newline at end of file