From ab270ebc9a2f59f7d232819d16ea021d5ccab8c9 Mon Sep 17 00:00:00 2001 From: randogoth Date: Tue, 4 Mar 2025 14:54:08 +0000 Subject: [PATCH] linstring approach. still fragmented --- water.py | 227 ++++++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 209 insertions(+), 18 deletions(-) diff --git a/water.py b/water.py index 5b97d3b..027eb5e 100644 --- a/water.py +++ b/water.py @@ -3,9 +3,143 @@ import argparse import pandas as pd import geopandas as gpd import geohash +import numpy as np import logging -from shapely.geometry import box -from shapely.ops import unary_union +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}") # Configure logging logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.INFO) @@ -33,8 +167,48 @@ def ensure_crs_consistency(gdf): 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 + + # 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"]] + + 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: removes old coastlines and inserts updated ones.""" + """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}...") @@ -55,29 +229,46 @@ def process_geohash_files(coastline_gdf, geohash_files): bbox = geohash_bbox(geohash_code) new_coastline = coastline_gdf[coastline_gdf.intersects(bbox)].copy() - # Clip coastline to the geohash boundary - new_coastline["geometry"] = new_coastline["geometry"].apply(lambda geom: geom.intersection(bbox)) - if new_coastline.empty: - logger.info(f"No coastline found for {geohash_code}. Skipping update.") + logger.info(f"No coastline found for {geohash_code}. Skipping coastline update.") continue - # Assign "natural: coastline" tag - new_coastline["tags"] = [{"natural": "coastline"}] * len(new_coastline) + # 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 - # Ensure columns match before merging - for col in ["feature_id", "tags"]: - if col not in filtered_gdf.columns: - filtered_gdf[col] = None - if col not in new_coastline.columns: - new_coastline[col] = None + if merged_coastline is None: + logger.warning(f"Geohash {geohash_code}: Coastline merging failed.") + continue - # Merge updated coastline into geohash file - updated_gdf = gpd.GeoDataFrame(pd.concat([filtered_gdf, new_coastline], ignore_index=True), crs="EPSG:4326") + logger.info(f"Adding water polygon for {geohash_code}...") + + # 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") # Save back to parquet updated_gdf.to_parquet(file_path, index=False) - logger.info(f"Updated coastline in {file_path}") + logger.info(f"Updated geohash file {file_path}.") def main(): parser = argparse.ArgumentParser(description="Update coastline geometries in existing geohash parquet files.")