feat(pygal): implement hexbin-map-geographic (#7746)

## Implementation: `hexbin-map-geographic` - python/pygal

Implements the **python/pygal** version of `hexbin-map-geographic`.

**File:** `plots/hexbin-map-geographic/implementations/python/pygal.py`

**Parent Issue:** #3767

---
:robot: *[impl-generate
workflow](https://github.com/MarkusNeusinger/anyplot/actions/runs/26513744161)*

---------

Co-authored-by: github-actions[bot] <github-actions[bot]@users.noreply.github.com>
Co-authored-by: Claude Sonnet 4.6 <noreply@anthropic.com>
Co-authored-by: Markus Neusinger <2921697+MarkusNeusinger@users.noreply.github.com>

Author: 3 people

plots/hexbin-map-geographic/implementations/python/pygal.py

@@ -1,18 +1,22 @@
-""" pyplots.ai
+""" anyplot.ai
 hexbin-map-geographic: Hexagonal Binning Map
-Library: pygal 3.1.0 | Python 3.13.11
-Quality: 75/100 | Created: 2026-01-20
+Library: pygal 3.1.0 | Python 3.13.13
+Quality: 84/100 | Updated: 2026-05-27
 """
 
-# Fix module name conflict (this file is named pygal.py)
+import math
+import os
+import re
 import sys
 from collections import defaultdict
 
 
+# Fix module name conflict (this file is named pygal.py)
 _cwd = sys.path[0] if sys.path and sys.path[0] else None
 if _cwd:
     sys.path.remove(_cwd)
 
+import cairosvg  # noqa: E402
 import numpy as np  # noqa: E402
 import pygal  # noqa: E402
 from pygal.style import Style  # noqa: E402
@@ -21,88 +25,85 @@
 if _cwd:
     sys.path.insert(0, _cwd)
 
-# Data - Simulated NYC taxi pickup locations (Manhattan area)
+# Theme tokens
+THEME = os.getenv("ANYPLOT_THEME", "light")
+PAGE_BG = "#FAF8F1" if THEME == "light" else "#1A1A17"
+INK = "#1A1A17" if THEME == "light" else "#F0EFE8"
+INK_MUTED = "#6B6A63" if THEME == "light" else "#A8A79F"
+
+# imprint_seq: #009E73 → #4467A3 (5 evenly-spaced stops, low → high density)
+seq_stops = ("#009E73", "#11907F", "#22838B", "#337597", "#4467A3")
+n_density_bins = 5
+
+# Data — NYC taxi pickup locations (Manhattan)
 np.random.seed(42)
 n_points = 5000
-
-# NYC Manhattan bounds (approximately)
 lat_min, lat_max = 40.70, 40.82
 lon_min, lon_max = -74.02, -73.93
 
-# Create multiple hotspots (Midtown, Lower Manhattan, Upper East Side)
-# Cluster 1: Midtown (Times Square area)
 c1_lat = np.random.normal(40.758, 0.015, n_points // 3)
 c1_lon = np.random.normal(-73.985, 0.01, n_points // 3)
-c1_vals = np.random.exponential(25, n_points // 3)  # Fare amounts
+c1_vals = np.random.exponential(25, n_points // 3)
 
-# Cluster 2: Lower Manhattan (Financial District)
 c2_lat = np.random.normal(40.710, 0.012, n_points // 3)
 c2_lon = np.random.normal(-74.010, 0.008, n_points // 3)
-c2_vals = np.random.exponential(35, n_points // 3)  # Higher fares downtown
+c2_vals = np.random.exponential(35, n_points // 3)
 
-# Cluster 3: Upper East Side
 c3_lat = np.random.normal(40.775, 0.018, n_points // 3)
 c3_lon = np.random.normal(-73.960, 0.012, n_points // 3)
 c3_vals = np.random.exponential(20, n_points // 3)
 
-# Combine all points
-lat = np.concatenate([c1_lat, c2_lat, c3_lat])
-lon = np.concatenate([c1_lon, c2_lon, c3_lon])
+lat = np.clip(np.concatenate([c1_lat, c2_lat, c3_lat]), lat_min, lat_max)
+lon = np.clip(np.concatenate([c1_lon, c2_lon, c3_lon]), lon_min, lon_max)
 values = np.concatenate([c1_vals, c2_vals, c3_vals])
 
-# Clip to bounds
-lat = np.clip(lat, lat_min, lat_max)
-lon = np.clip(lon, lon_min, lon_max)
-
-# Compute hexagonal binning with aggregation (inline - KISS principle)
+# Hexagonal binning with count and fare aggregation
 gridsize = 25
-x = np.asarray(lon)
-y = np.asarray(lat)
-
-# Compute data bounds
-x_min, x_max = x.min(), x.max()
-y_min = y.min()
-
-# Hexagon dimensions
-x_range = x_max - x_min
-hex_width = x_range / gridsize
+x_arr = np.asarray(lon)
+y_arr = np.asarray(lat)
+x_min_v, x_max_v = x_arr.min(), x_arr.max()
+y_min_v = y_arr.min()
+hex_width = (x_max_v - x_min_v) / gridsize
 hex_height = hex_width * np.sqrt(3) / 2
 
-# Convert points to hex grid coordinates with value aggregation
-bins = defaultdict(lambda: {"count": 0, "sum": 0.0, "values": []})
-
-for xi, yi, vi in zip(x, y, values, strict=True):
-    col = (xi - x_min) / hex_width
+bins = defaultdict(lambda: {"count": 0, "sum": 0.0})
+for xi, yi, vi in zip(x_arr, y_arr, values, strict=True):
+    col = (xi - x_min_v) / hex_width
     row_offset = (int(col) % 2) * 0.5
-    row = (yi - y_min) / hex_height - row_offset
+    row = (yi - y_min_v) / hex_height - row_offset
     col_idx = int(round(col))
     row_idx = int(round(row))
     bins[(col_idx, row_idx)]["count"] += 1
     bins[(col_idx, row_idx)]["sum"] += vi
-    bins[(col_idx, row_idx)]["values"].append(vi)
 
-# Convert bin indices back to coordinates with full statistics
 hex_data = []
-
 for (col_idx, row_idx), data in bins.items():
-    cx = x_min + col_idx * hex_width
+    cx = x_min_v + col_idx * hex_width
     row_offset = (col_idx % 2) * 0.5
-    cy = y_min + (row_idx + row_offset) * hex_height
+    cy = y_min_v + (row_idx + row_offset) * hex_height
     count = data["count"]
     total = data["sum"]
     mean = total / count if count > 0 else 0
     hex_data.append({"lon": cx, "lat": cy, "count": count, "sum": total, "mean": mean})
 
-# Extract arrays for plotting
-hex_lon = np.array([h["lon"] for h in hex_data])
-hex_lat = np.array([h["lat"] for h in hex_data])
 counts = np.array([h["count"] for h in hex_data])
 
-# Get count statistics for binning - use percentile-based bins for better distribution
-count_min, count_max = counts.min(), counts.max()
-count_range = count_max - count_min if count_max > count_min else 1
+# Percentile-based bin edges for balanced distribution
+bin_edges = np.percentile(counts, [0, 20, 40, 60, 80, 100])
+for i in range(1, len(bin_edges)):
+    if bin_edges[i] <= bin_edges[i - 1]:
+        bin_edges[i] = bin_edges[i - 1] + 1
+
+# Shortened labels to prevent legend truncation
+bin_labels = [
+    f"Low ({int(bin_edges[0])}–{int(bin_edges[1])})",
+    f"Med-Low ({int(bin_edges[1])}–{int(bin_edges[2])})",
+    f"Medium ({int(bin_edges[2])}–{int(bin_edges[3])})",
+    f"Med-High ({int(bin_edges[3])}–{int(bin_edges[4])})",
+    f"High ({int(bin_edges[4])}+)",
+]
 
-# Simplified Manhattan coastline (approximate outline)
+# Geographic outlines (Manhattan island + waterways)
 manhattan_outline = [
     (-74.020, 40.700),
     (-74.010, 40.705),
@@ -121,8 +122,6 @@
     (-74.015, 40.720),
     (-74.020, 40.700),
 ]
-
-# Hudson River (west boundary approximation)
 hudson_river = [
     (-74.035, 40.690),
     (-74.025, 40.705),
@@ -132,8 +131,6 @@
     (-73.985, 40.810),
     (-73.970, 40.835),
 ]
-
-# East River (east boundary approximation)
 east_river = [
     (-73.935, 40.695),
     (-73.940, 40.720),
@@ -143,44 +140,31 @@
     (-73.920, 40.835),
 ]
 
-# Custom style - YlOrRd colormap for density (5 levels)
-# Boundary lines use darker colors for visibility
+title = "hexbin-map-geographic · python · pygal · anyplot.ai"
+
 custom_style = Style(
-    background="white",
-    plot_background="#E8F4F8",  # Light water blue
-    foreground="#333333",
-    foreground_strong="#111111",
-    foreground_subtle="#666666",
-    guide_stroke_color="transparent",  # Hide grid lines completely
-    colors=(
-        # 3 boundary/river colors (darker for visibility)
-        "#555555",
-        "#446688",
-        "#446688",
-        # 5 density levels (YlOrRd colormap - darker for PNG visibility)
-        "#FFFFB2",  # Very low - light yellow
-        "#FECC5C",  # Low - yellow
-        "#FD8D3C",  # Medium - orange
-        "#F03B20",  # High - red-orange
-        "#BD0026",  # Very high - dark red
-    ),
-    opacity=0.85,  # Higher opacity for better visibility
+    background=PAGE_BG,
+    plot_background=PAGE_BG,
+    foreground=INK,
+    foreground_strong=INK,
+    foreground_subtle=INK_MUTED,
+    guide_stroke_color="transparent",
+    colors=(INK_MUTED, INK_MUTED, INK_MUTED) + seq_stops,
+    opacity=0.85,
     opacity_hover=0.95,
-    title_font_size=72,
-    label_font_size=48,
-    major_label_font_size=42,
-    legend_font_size=40,
+    title_font_size=66,
+    label_font_size=56,
+    major_label_font_size=44,
+    legend_font_size=44,
     value_font_size=36,
-    tooltip_font_size=36,
-    stroke_width=3,  # Thicker lines for boundaries
+    stroke_width=2.5,
 )
 
-# Create XY chart - disable grid for clean map background
 chart = pygal.XY(
-    width=4800,
-    height=2700,
+    width=3200,
+    height=1800,
     style=custom_style,
-    title="hexbin-map-geographic · pygal · pyplots.ai",
+    title=title,
     x_title="Longitude (°)",
     y_title="Latitude (°)",
     show_legend=True,
@@ -189,93 +173,74 @@
     legend_box_size=28,
     stroke=False,
     dots_size=30,
-    show_x_guides=False,  # Hide x grid lines
-    show_y_guides=False,  # Hide y grid lines
+    show_x_guides=False,
+    show_y_guides=False,
     explicit_size=True,
     print_values=False,
-    xrange=(lon_min - 0.015, lon_max + 0.015),  # Tighter margins
+    xrange=(lon_min - 0.015, lon_max + 0.015),
     range=(lat_min - 0.008, lat_max + 0.008),
 )
 
-# Add geographic boundaries as background lines with thicker stroke
+# Geographic boundaries (excluded from legend via None label)
 chart.add(None, manhattan_outline, stroke=True, dots_size=0, show_dots=False, fill=False, stroke_width=4)
 chart.add(None, hudson_river, stroke=True, dots_size=0, show_dots=False, fill=False, stroke_width=3)
 chart.add(None, east_river, stroke=True, dots_size=0, show_dots=False, fill=False, stroke_width=3)
 
-# Use percentile-based bin edges to ensure all 5 bins have data
-n_bins = 5
-percentiles = [0, 20, 40, 60, 80, 100]
-bin_edges = np.percentile(counts, percentiles)
-# Ensure edges are strictly increasing
-for i in range(1, len(bin_edges)):
-    if bin_edges[i] <= bin_edges[i - 1]:
-        bin_edges[i] = bin_edges[i - 1] + 1
-
-# Create legend labels showing density scale
-bin_labels = [
-    f"Low ({int(bin_edges[0])}-{int(bin_edges[1])})",
-    f"Medium-Low ({int(bin_edges[1])}-{int(bin_edges[2])})",
-    f"Medium ({int(bin_edges[2])}-{int(bin_edges[3])})",
-    f"Medium-High ({int(bin_edges[3])}-{int(bin_edges[4])})",
-    f"High ({int(bin_edges[4])}+)",
-]
-
-# Create series for each density level with rich tooltips
-series_data = [[] for _ in range(n_bins)]
+# Density bins with size encoding (larger hexagons = more pickups)
+series_data = [[] for _ in range(n_density_bins)]
 for h in hex_data:
     hx, hy = h["lon"], h["lat"]
     count = h["count"]
     total = h["sum"]
     mean = h["mean"]
-    # Assign to bin based on percentile edges
     bin_idx = 0
-    for i in range(1, n_bins):
+    for i in range(1, n_density_bins):
         if count >= bin_edges[i]:
             bin_idx = i
-    # Rich tooltip with cell statistics as spec requires
-    tooltip = (
-        f"Count: {count} pickups | Fares: ${total:.0f} total, ${mean:.2f} avg | Coords: ({hy:.4f}°N, {abs(hx):.4f}°W)"
-    )
-    point = {"value": (float(hx), float(hy)), "label": tooltip}
-    series_data[bin_idx].append(point)
-
-# Add ALL series (even empty ones) to ensure complete legend
-# Using larger dot sizes for density visualization
+    tooltip = f"Count: {count} | Fares: ${total:.0f} total, ${mean:.2f} avg | ({hy:.4f}°N, {abs(hx):.4f}°W)"
+    series_data[bin_idx].append({"value": (float(hx), float(hy)), "label": tooltip})
+
 dot_sizes = [26, 34, 44, 54, 66]
-for i in range(n_bins):
-    # Add dummy point off-screen if series is empty to show in legend
+for i in range(n_density_bins):
     if not series_data[i]:
-        # Add invisible point outside plot range to ensure legend appears
         series_data[i].append({"value": (-99, 0), "label": "No data"})
     chart.add(bin_labels[i], series_data[i], dots_size=dot_sizes[i])
 
-# Save PNG (primary output)
-chart.render_to_png("plot.png")
-
-# Save interactive HTML with CSS hexagon styling
-# Use CSS clip-path to transform circles into hexagons in interactive view
-html_content = f"""<!DOCTYPE html>
-<html>
-<head>
-    <meta charset="utf-8">
-    <title>hexbin-map-geographic - pygal</title>
-    <style>
-        body {{ margin: 0; display: flex; justify-content: center; align-items: center;
-               min-height: 100vh; background: #f5f5f5; }}
-        .chart {{ max-width: 100%; height: auto; }}
-        /* Transform dots into hexagons using CSS clip-path */
-        .dot {{
-            clip-path: polygon(50% 0%, 100% 25%, 100% 75%, 50% 100%, 0% 75%, 0% 25%);
-        }}
-    </style>
-</head>
-<body>
-    <figure class="chart">
-        {chart.render(is_unicode=True)}
-    </figure>
-</body>
-</html>
-"""
 
-with open("plot.html", "w", encoding="utf-8") as f:
-    f.write(html_content)
+def circles_to_hexagons(svg_text):
+    """Post-process SVG: replace circular dot markers with flat-top hexagonal polygons."""
+
+    def replace_one(m):
+        tag = m.group(0)
+        cx_m = re.search(r'\bcx="([^"]+)"', tag)
+        cy_m = re.search(r'\bcy="([^"]+)"', tag)
+        r_m = re.search(r'\br="([^"]+)"', tag)
+        if not (cx_m and cy_m and r_m):
+            return tag
+        cx = float(cx_m.group(1))
+        cy = float(cy_m.group(1))
+        r = float(r_m.group(1))
+        # Flat-top hexagon: vertex angles 30°, 90°, 150°, 210°, 270°, 330°
+        pts = " ".join(
+            f"{cx + r * math.cos(math.radians(60 * k + 30)):.2f},{cy + r * math.sin(math.radians(60 * k + 30)):.2f}"
+            for k in range(6)
+        )
+        poly = tag.replace("<circle", "<polygon", 1)
+        poly = re.sub(r'\bcx="[^"]*"\s*', "", poly)
+        poly = re.sub(r'\bcy="[^"]*"\s*', "", poly)
+        poly = re.sub(r'\br="[^"]*"', f'points="{pts}"', poly)
+        return poly
+
+    return re.sub(r"<circle\b[^>]*/>", replace_one, svg_text)
+
+
+# Render SVG, convert circle markers to hexagonal polygons, then produce PNG
+svg_bytes = chart.render()
+svg_str = svg_bytes.decode("utf-8")
+modified_svg = circles_to_hexagons(svg_str)
+modified_bytes = modified_svg.encode("utf-8")
+
+cairosvg.svg2png(bytestring=modified_bytes, write_to=f"plot-{THEME}.png")
+
+with open(f"plot-{THEME}.html", "wb") as f:
+    f.write(modified_bytes)