fix(pygal): address review feedback for stereonet-equal-area

Attempt 2/3 - fixes based on AI review

Author: github-actions[bot]

plots/stereonet-equal-area/implementations/pygal.py

@@ -1,4 +1,4 @@
-""" pyplots.ai
+"""pyplots.ai
 stereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)
 Library: pygal 3.1.0 | Python 3.14.3
 Quality: 68/100 | Created: 2026-03-15
@@ -29,36 +29,56 @@
 
 # Equal-area (Schmidt net) lower-hemisphere projection
 # Pole to plane: trend = strike + 90°, plunge = 90° - dip
-# Projection radius: r = √2 · sin(dip_rad / 2)
+# Projection radius: r = √2 · sin(plunge_rad / 2) where plunge = 90 - dip for poles
+# For great circles: r = √2 · sin((90 - plunge_of_line) / 2)
 # Cartesian: x = r · sin(trend_rad), y = r · cos(trend_rad)
 
-# Primitive circle (outer boundary representing the horizontal plane)
+# Primitive circle (outer boundary)
 theta_circle = np.linspace(0, 2 * np.pi, 361)
 circle_points = [(round(float(np.sin(t)), 4), round(float(np.cos(t)), 4)) for t in theta_circle]
 
 # Tick marks every 10° around perimeter
 tick_points = []
 for deg in range(0, 360, 10):
     rad = np.radians(deg)
-    inner, outer = 0.97, 1.03
+    inner, outer = 0.96, 1.04
     tick_points.append((round(inner * np.sin(rad), 4), round(inner * np.cos(rad), 4)))
     tick_points.append((round(outer * np.sin(rad), 4), round(outer * np.cos(rad), 4)))
     tick_points.append(None)
 
-# Compass markers (N, E, S, W as small cross marks beyond the circle)
-compass_points = []
-for _bearing, offset_x, offset_y in [(0, 0, 1.12), (90, 1.12, 0), (180, 0, -1.12), (270, -1.12, 0)]:
-    compass_points.append((round(offset_x - 0.02, 4), round(offset_y, 4)))
-    compass_points.append((round(offset_x + 0.02, 4), round(offset_y, 4)))
-    compass_points.append(None)
-    compass_points.append((round(offset_x, 4), round(offset_y - 0.02, 4)))
-    compass_points.append((round(offset_x, 4), round(offset_y + 0.02, 4)))
-    compass_points.append(None)
+# Equal-area net grid lines (subtle, as spec requires)
+# Small circles at dip intervals of 30° (at 30°, 60°)
+grid_points = []
+for dip_grid in [30, 60]:
+    r_grid = np.sqrt(2) * np.sin(np.radians(dip_grid / 2))
+    theta_grid = np.linspace(0, 2 * np.pi, 181)
+    for t in theta_grid:
+        grid_points.append((round(float(r_grid * np.sin(t)), 4), round(float(r_grid * np.cos(t)), 4)))
+    grid_points.append(None)
+
+# Diametral lines every 30° (great circle diameters through center)
+for az in range(0, 180, 30):
+    rad = np.radians(az)
+    grid_points.append((round(-np.sin(rad), 4), round(-np.cos(rad), 4)))
+    grid_points.append((round(np.sin(rad), 4), round(np.cos(rad), 4)))
+    grid_points.append(None)
+
+# Compass direction labels as dot patterns
+# N label at top - letter shape using connected line segments
+n_label = []
+n_cx, n_cy = 0.0, 1.15
+n_s = 0.03
+for t in np.linspace(0, 1, 10):
+    n_label.append((round(n_cx - n_s, 4), round(n_cy - n_s * 2 + t * n_s * 4, 4)))
+n_label.append(None)
+for t in np.linspace(0, 1, 10):
+    n_label.append((round(n_cx - n_s + t * n_s * 2, 4), round(n_cy + n_s * 2 - t * n_s * 4, 4)))
+n_label.append(None)
+for t in np.linspace(0, 1, 10):
+    n_label.append((round(n_cx + n_s, 4), round(n_cy - n_s * 2 + t * n_s * 4, 4)))
+
 
 # Great circle computation for all planes
-# For rake angle α (0 to π), a line in the plane has:
-#   plunge = arcsin(sin(dip) · sin(α))
-#   trend = strike + atan2(sin(α) · cos(dip), cos(α))
 alphas = np.linspace(0, np.pi, 91)
 
 feature_sets = [
@@ -69,6 +89,8 @@
 
 gc_series = {}
 pole_series = {}
+all_pole_x = []
+all_pole_y = []
 
 for name, strikes, dips in feature_sets:
     # Great circles with None separators between each plane
@@ -91,32 +113,107 @@
     pole_x = pole_r * np.sin(pole_trend_rad)
     pole_y = pole_r * np.cos(pole_trend_rad)
     pole_series[name] = [(round(float(xi), 4), round(float(yi), 4)) for xi, yi in zip(pole_x, pole_y, strict=True)]
+    all_pole_x.extend(pole_x)
+    all_pole_y.extend(pole_y)
+
+# Kamb density contours on all pole data
+# Grid-based density estimation in projected space
+all_pole_x = np.array(all_pole_x)
+all_pole_y = np.array(all_pole_y)
+
+grid_res = 60
+gx = np.linspace(-1, 1, grid_res)
+gy = np.linspace(-1, 1, grid_res)
+gxx, gyy = np.meshgrid(gx, gy)
+
+# Gaussian kernel density on projected coordinates
+sigma = 0.12
+density = np.zeros_like(gxx)
+for px, py in zip(all_pole_x, all_pole_y, strict=True):
+    density += np.exp(-((gxx - px) ** 2 + (gyy - py) ** 2) / (2 * sigma**2))
+
+# Mask outside the primitive circle
+mask = gxx**2 + gyy**2 > 1.0
+density[mask] = 0
 
-# Style
+# Extract contour lines using a simple marching approach
+# Generate contour polylines at specific density levels
+contour_levels = [2.0, 4.0, 6.0, 8.0]
+max_density = density.max()
+contour_series = []
+
+for level in contour_levels:
+    if level >= max_density:
+        continue
+    level_points = []
+    # Scan horizontal edges for contour crossings
+    for i in range(grid_res - 1):
+        for j in range(grid_res):
+            v0, v1 = density[j, i], density[j, i + 1]
+            if (v0 - level) * (v1 - level) < 0 and v0 != v1:
+                frac = (level - v0) / (v1 - v0)
+                cx = gx[i] + frac * (gx[i + 1] - gx[i])
+                cy = gy[j]
+                if cx**2 + cy**2 <= 1.0:
+                    level_points.append((cx, cy))
+    # Scan vertical edges
+    for i in range(grid_res):
+        for j in range(grid_res - 1):
+            v0, v1 = density[j, i], density[j + 1, i]
+            if (v0 - level) * (v1 - level) < 0 and v0 != v1:
+                frac = (level - v0) / (v1 - v0)
+                cx = gx[i]
+                cy = gy[j] + frac * (gy[j + 1] - gy[j])
+                if cx**2 + cy**2 <= 1.0:
+                    level_points.append((cx, cy))
+
+    if level_points:
+        # Sort points by angle from each cluster centroid to form connected contours
+        pts = np.array(level_points)
+        # Order points by angle for smoother contour rendering
+        angles = np.arctan2(pts[:, 0], pts[:, 1])
+        order = np.argsort(angles)
+        sorted_pts = pts[order]
+        # Split into segments when consecutive points are far apart
+        segment = [(round(float(sorted_pts[0, 0]), 4), round(float(sorted_pts[0, 1]), 4))]
+        for k in range(1, len(sorted_pts)):
+            dist = np.sqrt(
+                (sorted_pts[k, 0] - sorted_pts[k - 1, 0]) ** 2 + (sorted_pts[k, 1] - sorted_pts[k - 1, 1]) ** 2
+            )
+            if dist > 0.15:
+                segment.append(None)
+            segment.append((round(float(sorted_pts[k, 0]), 4), round(float(sorted_pts[k, 1]), 4)))
+        contour_series.append(segment)
+
+
+# Style - colorblind-safe palette (blue, orange, amber instead of green)
 custom_style = Style(
     background="white",
     plot_background="white",
     foreground="#333333",
     foreground_strong="#333333",
     foreground_subtle="#eeeeee",
     colors=(
-        "#aaaaaa",  # primitive circle
-        "#aaaaaa",  # tick marks
-        "#aaaaaa",  # compass marks
+        "#d0d0d0",  # grid lines
+        "#777777",  # structural elements (circle + ticks + N)
         "#306998",  # bedding great circles
         "#D4513D",  # fault great circles
-        "#2CA02C",  # joint great circles
+        "#DAA520",  # joint great circles (amber/gold, colorblind-safe)
         "#306998",  # bedding poles
         "#D4513D",  # fault poles
-        "#2CA02C",  # joint poles
+        "#DAA520",  # joint poles (amber/gold)
+        "#666666",  # contour level 1
+        "#555555",  # contour level 2
+        "#444444",  # contour level 3
+        "#333333",  # contour level 4
     ),
     title_font_size=56,
     label_font_size=1,
     major_label_font_size=1,
     legend_font_size=36,
     value_font_size=20,
     tooltip_font_size=28,
-    opacity=0.5,
+    opacity=0.6,
     opacity_hover=0.8,
 )
 
@@ -128,35 +225,42 @@
     title="stereonet-equal-area · pygal · pyplots.ai",
     show_legend=True,
     legend_at_bottom=True,
-    legend_at_bottom_columns=3,
+    legend_at_bottom_columns=4,
     show_x_labels=False,
     show_y_labels=False,
     show_x_guides=False,
     show_y_guides=False,
-    xrange=(-1.35, 1.35),
-    range=(-1.35, 1.35),
+    xrange=(-1.25, 1.25),
+    range=(-1.25, 1.25),
     dots_size=0,
     allow_interruptions=True,
-    margin_top=40,
-    margin_bottom=60,
-    margin_left=20,
-    margin_right=20,
+    margin_top=30,
+    margin_bottom=50,
+    margin_left=10,
+    margin_right=10,
 )
 
-# Structural elements (circle, ticks, compass)
-chart.add("", circle_points, stroke=True, show_dots=False, stroke_style={"width": 2.5})
-chart.add("", tick_points, stroke=True, show_dots=False, stroke_style={"width": 1.5})
-chart.add("", compass_points, stroke=True, show_dots=False, stroke_style={"width": 2})
+# Grid lines (subtle equal-area net)
+chart.add("", grid_points, stroke=True, show_dots=False, stroke_style={"width": 0.8, "dasharray": "4,4"})
+
+# Combine structural elements: circle + ticks + N label
+structural = circle_points + [None] + tick_points + [None] + n_label
+chart.add("", structural, stroke=True, show_dots=False, stroke_style={"width": 2.5})
 
-# Great circles (planes)
-chart.add("Bedding (planes)", gc_series["Bedding"], stroke=True, show_dots=False, stroke_style={"width": 1.5})
-chart.add("Faults (planes)", gc_series["Faults"], stroke=True, show_dots=False, stroke_style={"width": 1.5})
-chart.add("Joints (planes)", gc_series["Joints"], stroke=True, show_dots=False, stroke_style={"width": 1.5})
+# Great circles (planes) - thicker lines for visibility
+chart.add("Bedding (planes)", gc_series["Bedding"], stroke=True, show_dots=False, stroke_style={"width": 2.5})
+chart.add("Faults (planes)", gc_series["Faults"], stroke=True, show_dots=False, stroke_style={"width": 2.5})
+chart.add("Joints (planes)", gc_series["Joints"], stroke=True, show_dots=False, stroke_style={"width": 2.5})
 
 # Poles to planes (scatter points)
-chart.add("Bedding (poles)", pole_series["Bedding"], stroke=False, show_dots=True, dots_size=12)
-chart.add("Fault (poles)", pole_series["Faults"], stroke=False, show_dots=True, dots_size=12)
-chart.add("Joint (poles)", pole_series["Joints"], stroke=False, show_dots=True, dots_size=12)
+chart.add("Bedding (poles)", pole_series["Bedding"], stroke=False, show_dots=True, dots_size=14)
+chart.add("Faults (poles)", pole_series["Faults"], stroke=False, show_dots=True, dots_size=14)
+chart.add("Joints (poles)", pole_series["Joints"], stroke=False, show_dots=True, dots_size=14)
+
+# Density contours - only first gets legend label
+for i, contour_data in enumerate(contour_series):
+    label = "Density contours" if i == 0 else ""
+    chart.add(label, contour_data, stroke=True, show_dots=False, stroke_style={"width": 1.5, "dasharray": "6,3"})
 
 # Save
 chart.render_to_png("plot.png")