diff --git a/create_map_poster.py b/create_map_poster.py index 3eab412b..b7aa6d04 100755 --- a/create_map_poster.py +++ b/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) diff --git a/posters/durban_terracotta_20260202_223000.png b/posters/durban_terracotta_20260202_223000.png new file mode 100644 index 00000000..b704cf66 Binary files /dev/null and b/posters/durban_terracotta_20260202_223000.png differ diff --git a/posters/hong_kong_terracotta_20260210_211635.png b/posters/hong_kong_terracotta_20260210_211635.png new file mode 100644 index 00000000..c9bc658a Binary files /dev/null and b/posters/hong_kong_terracotta_20260210_211635.png differ diff --git a/posters/san_francisco_terracotta_20260210_211551.png b/posters/san_francisco_terracotta_20260210_211551.png new file mode 100644 index 00000000..abc8529a Binary files /dev/null and b/posters/san_francisco_terracotta_20260210_211551.png differ