Feature/renderoceans

create_map_poster.py

@@ -22,12 +22,13 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import osmnx as ox
-from geopandas import GeoDataFrame
+from geopandas import GeoDataFrame, GeoSeries
 from geopy.geocoders import Nominatim
 from lat_lon_parser import parse
 from matplotlib.font_manager import FontProperties
 from networkx import MultiDiGraph
-from shapely.geometry import Point
+from shapely.geometry import Point, box
+from shapely.ops import linemerge, polygonize, unary_union
 from tqdm import tqdm
 
 from font_management import load_fonts
@@ -493,6 +494,7 @@ def create_poster(
     display_city=None,
     display_country=None,
     fonts=None,
+    include_oceans=True,
 ):
     """
     Generate a complete map poster with roads, water, parks, and typography.
@@ -542,7 +544,12 @@ def create_poster(
         water = fetch_features(
             point,
             compensated_dist,
-            tags={"natural": ["water", "bay", "strait"], "waterway": "riverbank"},
+            tags={
+                "natural": ["water", "bay", "strait", "coastline", "sea", "ocean"],
+                "waterway": ["riverbank", "river", "canal", "dock", "basin"],
+                "place": ["sea", "ocean", "bay"],
+                "water": True
+            },
             name="water",
         )
         pbar.update(1)
@@ -568,36 +575,117 @@ def create_poster(
     # Project graph to a metric CRS so distances and aspect are linear (meters)
     g_proj = ox.project_graph(g)
 
+    # Determine cropping limits to maintain the poster aspect ratio
+    # We do this early so we can use it for coastline polygonization
+    crop_xlim, crop_ylim = get_crop_limits(g_proj, point, fig, compensated_dist)
+
     # 3. Plot Layers
-    # Layer 1: Polygons (filter to only plot polygon/multipolygon geometries, not points)
+    # Layer 1: Water
     if water is not None and not water.empty:
-        # Filter to only polygon/multipolygon geometries to avoid point features showing as dots
-        water_polys = water[water.geometry.type.isin(["Polygon", "MultiPolygon"])]
-        if not water_polys.empty:
-            # Project water features in the same CRS as the graph
-            try:
-                water_polys = ox.projection.project_gdf(water_polys)
-            except Exception:
-                water_polys = water_polys.to_crs(g_proj.graph['crs'])
-            water_polys.plot(ax=ax, facecolor=THEME['water'], edgecolor='none', zorder=0.5)
+        # Project water features in the same CRS as the graph
+        water_proj = water.to_crs(g_proj.graph['crs'])
+
+        # Separate water polygons and island polygons
+        water_polys_all = water_proj[water_proj.geometry.type.isin(["Polygon", "MultiPolygon"])]
+        if not water_polys_all.empty:
+            # Filter out islands/islets that might be tagged with natural=coastline, These should be treated as land.
+            is_island = np.array([False] * len(water_polys_all))
+            if "place" in water_polys_all.columns:
+                is_island = np.array(water_polys_all["place"].isin(["island", "islet", "archipelago"]))
+
+            water_polys = water_polys_all[~is_island]
+            island_polys = water_polys_all[is_island]
+
+            if not water_polys.empty:
+                water_polys.plot(ax=ax, facecolor=THEME['water'], edgecolor='none', zorder=0.5)
+
+            if not island_polys.empty:
+                # Plot islands with background color to ensure they aren't covered by ocean fill
+                island_polys.plot(ax=ax, facecolor=THEME['bg'], edgecolor='none', zorder=0.6)
+
+        # Handle coastlines - try to form filled ocean polygons
+        water_lines = water_proj[water_proj.geometry.type.isin(["LineString", "MultiLineString"])]
+        if include_oceans and not water_lines.empty:
+            # Look for lines that represent coastlines or boundaries
+            coast_tags = ["coastline", "bay", "strait"]
+            if "natural" in water_lines.columns:
+                coastlines = water_lines[water_lines["natural"].isin(coast_tags)]
+
+                if not coastlines.empty:
+                    bbox_poly = box(crop_xlim[0], crop_ylim[0], crop_xlim[1], crop_ylim[1])
+                    merged_coast = linemerge(list(coastlines.geometry.values))
+                    if not merged_coast.is_empty:
+                        # Intersect with a slightly buffered bbox to ensure we catch edges
+                        merged_coast = merged_coast.intersection(bbox_poly.buffer(10))
+
+                        # Union with bbox boundary to form closed areas
+                        combined = unary_union([merged_coast, bbox_poly.boundary])
+                        candidate_polys = list(polygonize(combined))
+
+                        if candidate_polys:
+                            # Heuristic: Land polygons contain many road nodes, Water polygons contain few/none
+                            # Get road nodes as points
+                            node_points = [Point(d['x'], d['y']) for n, d in g_proj.nodes(data=True)]
+
+                            if node_points:
+                                sample_count = min(800, len(node_points))
+                                sample_nodes = node_points[::max(1, len(node_points)//sample_count)]
+
+                                # Also get park centers to use as "known land" points
+                                park_points = []
+                                if parks is not None and not parks.empty:
+                                    try:
+                                        parks_proj = ox.projection.project_gdf(parks)
+                                        park_points = [p.centroid for p in parks_proj.geometry if p is not None]
+                                    except Exception:
+                                        pass
+
+                                for poly in candidate_polys:
+                                    # Skip very small fragments
+                                    if poly.area < (bbox_poly.area * 0.001):
+                                        continue
+
+                                    # Count how many road nodes fall into this polygon
+                                    nodes_inside = sum(1 for p in sample_nodes if poly.contains(p))
+
+                                    # Heuristic: Land polygons have significantly higher road node density than water.
+                                    # We use density (nodes per km²) to be scale-invariant.
+                                    # A threshold of 50 nodes/km² (with sampling) is safe to distinguish 
+                                    # even sparsely populated land from water/oceans.
+
+                                    # Adjust density based on sampling rate
+                                    sampling_ratio = len(sample_nodes) / len(node_points)
+                                    estimated_total_nodes = nodes_inside / sampling_ratio if sampling_ratio > 0 else 0
+                                    density = (estimated_total_nodes / poly.area) * 1_000_000 # nodes per km²
+
+                                    is_water = False
+                                    if density < 50:
+                                        is_water = True
+
+                                    if is_water:
+                                        GeoSeries([poly]).plot(ax=ax, facecolor=THEME['water'], edgecolor='none', zorder=0.5)
+
+        if not water_lines.empty:
+            lines_to_plot = water_lines
+            # If oceans are disabled, filter out coastline outlines
+            if not include_oceans and "natural" in lines_to_plot.columns:
+                lines_to_plot = lines_to_plot[lines_to_plot["natural"] != "coastline"]
+
+            if not lines_to_plot.empty:
+                lines_to_plot.plot(ax=ax, color=THEME['water'], linewidth=0.8, zorder=0.5)
 
     if parks is not None and not parks.empty:
         # Filter to only polygon/multipolygon geometries to avoid point features showing as dots
         parks_polys = parks[parks.geometry.type.isin(["Polygon", "MultiPolygon"])]
         if not parks_polys.empty:
             # Project park features in the same CRS as the graph
-            try:
-                parks_polys = ox.projection.project_gdf(parks_polys)
-            except Exception:
-                parks_polys = parks_polys.to_crs(g_proj.graph['crs'])
+            parks_polys = parks_polys.to_crs(g_proj.graph['crs'])
             parks_polys.plot(ax=ax, facecolor=THEME['parks'], edgecolor='none', zorder=0.8)
     # Layer 2: Roads with hierarchy coloring
     print("Applying road hierarchy colors...")
     edge_colors = get_edge_colors_by_type(g_proj)
     edge_widths = get_edge_widths_by_type(g_proj)
 
-    # Determine cropping limits to maintain the poster aspect ratio
-    crop_xlim, crop_ylim = get_crop_limits(g_proj, point, fig, compensated_dist)
     # Plot the projected graph and then apply the cropped limits
     ox.plot_graph(
         g_proj, ax=ax, bgcolor=THEME['bg'],
@@ -955,6 +1043,13 @@ def list_themes():
         choices=["png", "svg", "pdf"],
         help="Output format for the poster (default: png)",
     )
+    parser.add_argument(
+        "--include-oceans",
+        dest="include_oceans",
+        action="store_true",
+        help="Enable automatic ocean/sea filling based on coastlines",
+    )
+    parser.set_defaults(include_oceans=False)
 
     args = parser.parse_args()
 
@@ -1037,6 +1132,7 @@ def list_themes():
                 display_city=args.display_city,
                 display_country=args.display_country,
                 fonts=custom_fonts,
+                include_oceans=args.include_oceans,
             )
 
         print("\n" + "=" * 50)