{"spec_id":"stereonet-equal-area","library":"pygal","language":"python","code":"\"\"\" anyplot.ai\nstereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\nLibrary: pygal 3.1.0 | Python 3.13.13\nQuality: 90/100 | Updated: 2026-06-16\n\"\"\"\n\nimport os\nimport sys\n\n\n# Remove this file's directory from sys.path so `import pygal` resolves to the\n# installed package, not this script (which shares the same name).\n_this_dir = os.path.dirname(os.path.abspath(__file__))\nsys.path = [p for p in sys.path if os.path.abspath(p or \".\") != _this_dir]\n\nimport numpy as np\nimport pygal\nfrom pygal.style import Style\n\n\n# Theme-adaptive chrome (Imprint palette) ------------------------------------\nTHEME = os.getenv(\"ANYPLOT_THEME\", \"light\")\nPAGE_BG = \"#FAF8F1\" if THEME == \"light\" else \"#1A1A17\"\nINK = \"#1A1A17\" if THEME == \"light\" else \"#F0EFE8\"\nINK_SOFT = \"#4A4A44\" if THEME == \"light\" else \"#B8B7B0\"\nINK_MUTED = \"#6B6A63\" if THEME == \"light\" else \"#A8A79F\"\n\n# Imprint categorical palette positions used here:\nBEDDING = \"#009E73\"  # brand green — ALWAYS first series\nFAULT = \"#C475FD\"  # lavender\nJOINT = \"#4467A3\"  # blue\nDENSITY = \"#AE3030\"  # matte red — pole-concentration \"hotspot\" emphasis\n\n# Data - field measurements from a geological mapping campaign ----------------\nnp.random.seed(42)\n\n# Bedding planes - consistent NE strike (~040 deg), moderate SE dip (~30 deg)\nn_bedding = 25\nbedding_strike = np.random.normal(40, 8, n_bedding) % 360\nbedding_dip = np.clip(np.random.normal(30, 5, n_bedding), 5, 85)\n\n# Fault planes - steeper, striking ESE (~120 deg), steep dip (~65 deg)\nn_faults = 15\nfault_strike = np.random.normal(120, 12, n_faults) % 360\nfault_dip = np.clip(np.random.normal(65, 8, n_faults), 10, 89)\n\n# Joint set - sub-vertical (~80 deg), striking roughly N-S (~350 deg)\nn_joints = 20\njoint_strike = np.random.normal(350, 10, n_joints) % 360\njoint_dip = np.clip(np.random.normal(80, 6, n_joints), 15, 89)\n\n# Equal-area (Schmidt net) lower-hemisphere projection.\n# Pole to plane: trend = strike + 90 deg, plunge = 90 deg - dip.\n# Equal-area radius for a line of plunge p: r = sqrt(2) * sin((90 - p) / 2).\n\n\ndef clip_to_circle(x, y):\n    \"\"\"Clip a point to the primitive circle (r <= 1).\"\"\"\n    r = np.hypot(x, y)\n    return (x / r, y / r) if r > 1.0 else (x, y)\n\n\n# Great circles + poles for each feature set ---------------------------------\nalphas = np.linspace(0, np.pi, 91)\nfeature_sets = [\n    (\"Bedding\", bedding_strike, bedding_dip),\n    (\"Faults\", fault_strike, fault_dip),\n    (\"Joints\", joint_strike, joint_dip),\n]\n\n# Visual hierarchy: bedding primary (thickest), faults secondary, joints tertiary\ngc_widths = {\"Bedding\": 3.0, \"Faults\": 2.2, \"Joints\": 1.6}\npole_sizes = {\"Bedding\": 16, \"Faults\": 12, \"Joints\": 9}\n\ngc_series, pole_series = {}, {}\nall_pole_x, all_pole_y = [], []\n\nfor name, strikes, dips in feature_sets:\n    gc_data = []\n    for s, d in zip(strikes, dips, strict=True):\n        d_rad = np.radians(d)\n        plunges = np.degrees(np.arcsin(np.sin(d_rad) * np.sin(alphas)))\n        trends = s + np.degrees(np.arctan2(np.sin(alphas) * np.cos(d_rad), np.cos(alphas)))\n        r = np.sqrt(2) * np.sin(np.radians((90 - plunges) / 2))\n        x = r * np.sin(np.radians(trends))\n        y = r * np.cos(np.radians(trends))\n        if gc_data:\n            gc_data.append(None)\n        for xi, yi in zip(x, y, strict=True):\n            cx, cy = clip_to_circle(float(xi), float(yi))\n            gc_data.append((round(cx, 4), round(cy, 4)))\n    gc_series[name] = gc_data\n\n    # Poles to planes (normal to each plane)\n    pole_trend = np.radians((strikes + 90) % 360)\n    pole_r = np.sqrt(2) * np.sin(np.radians(dips / 2))\n    px = pole_r * np.sin(pole_trend)\n    py = pole_r * np.cos(pole_trend)\n    pole_series[name] = [(round(float(a), 4), round(float(b), 4)) for a, b in zip(px, py, strict=True)]\n    all_pole_x.extend(px)\n    all_pole_y.extend(py)\n\n# Pole-density estimate + contour extraction ---------------------------------\n# Gaussian density on the projected poles. pygal has no contour primitive, so a\n# small marching-squares pass emits the iso-level crossing segments directly —\n# pygal's allow_interruptions then strings the disconnected segments together.\nall_pole_x = np.array(all_pole_x)\nall_pole_y = np.array(all_pole_y)\n\ngrid_res = 110\ngx = gy = np.linspace(-1, 1, grid_res)\ngxx, gyy = np.meshgrid(gx, gy)\nsigma = 0.12\ndensity = np.zeros_like(gxx)\nfor px, py in zip(all_pole_x, all_pole_y, strict=True):\n    density += np.exp(-((gxx - px) ** 2 + (gyy - py) ** 2) / (2 * sigma**2))\ndensity[gxx**2 + gyy**2 > 1.0] = 0.0  # mask outside the primitive circle\n\n\ndef contour_segments(d, level):\n    \"\"\"Marching-squares iso-line segments at `level` (no polyline chaining).\"\"\"\n    segs = []\n    corners = ((0, 1), (1, 2), (2, 3), (3, 0))  # cell edges as corner index pairs\n    for j in range(d.shape[0] - 1):\n        for i in range(d.shape[1] - 1):\n            v = (d[j, i], d[j, i + 1], d[j + 1, i + 1], d[j + 1, i])\n            above = tuple(val >= level for val in v)\n            if all(above) or not any(above):\n                continue\n            pos = ((gx[i], gy[j]), (gx[i + 1], gy[j]), (gx[i + 1], gy[j + 1]), (gx[i], gy[j + 1]))\n            cross = {}\n            for a, b in corners:\n                if above[a] != above[b]:\n                    t = 0.5 if abs(v[b] - v[a]) < 1e-12 else (level - v[a]) / (v[b] - v[a])\n                    cross[(a, b)] = (pos[a][0] + t * (pos[b][0] - pos[a][0]), pos[a][1] + t * (pos[b][1] - pos[a][1]))\n            pts = list(cross.values())\n            if len(pts) == 2:\n                segs.append((pts[0], pts[1]))\n            elif len(pts) == 4:  # saddle - pair edges by the cell-centre value\n                e = cross\n                pairs = [(0, 1), (2, 3)] if sum(v) / 4 >= level else [(0, 3), (1, 2)]\n                keys = list(e)\n                segs += [(e[keys[a]], e[keys[b]]) for a, b in pairs]\n    return segs\n\n\ncontour_data = []\nfor lvl in (lvl for lvl in (2.0, 4.0, 6.0, 8.0) if lvl < density.max()):\n    for p0, p1 in contour_segments(density, lvl):\n        if contour_data:\n            contour_data.append(None)\n        x0, y0 = clip_to_circle(*p0)\n        x1, y1 = clip_to_circle(*p1)\n        contour_data.append((round(x0, 4), round(y0, 4)))\n        contour_data.append((round(x1, 4), round(y1, 4)))\n\n# Equal-area net grid (subtle) -----------------------------------------------\ngrid_data = []\nfor dip_grid in (30, 60):  # small circles at 30 / 60 deg dip\n    rg = np.sqrt(2) * np.sin(np.radians(dip_grid / 2))\n    for t in np.linspace(0, 2 * np.pi, 181):\n        grid_data.append((round(float(rg * np.sin(t)), 4), round(float(rg * np.cos(t)), 4)))\n    grid_data.append(None)\nfor az in range(0, 180, 30):  # diametral lines every 30 deg\n    rad = np.radians(az)\n    grid_data.append((round(-np.sin(rad), 4), round(-np.cos(rad), 4)))\n    grid_data.append((round(np.sin(rad), 4), round(np.cos(rad), 4)))\n    grid_data.append(None)\n\n# Primitive circle + perimeter ticks + cardinal letters ----------------------\nboundary_data = []\nfor t in np.linspace(0, 2 * np.pi, 361):  # primitive circle\n    boundary_data.append((round(float(np.sin(t)), 4), round(float(np.cos(t)), 4)))\nboundary_data.append(None)\nfor deg in range(0, 360, 10):  # ticks every 10 deg\n    rad = np.radians(deg)\n    inner, outer = (0.97, 1.05) if deg % 30 == 0 else (0.99, 1.03)\n    boundary_data.append((round(inner * np.sin(rad), 4), round(inner * np.cos(rad), 4)))\n    boundary_data.append((round(outer * np.sin(rad), 4), round(outer * np.cos(rad), 4)))\n    boundary_data.append(None)\n\n# Cardinal labels (N/E/S/W) drawn as polyline letters just outside the circle\nLETTER_STROKES = {\n    \"N\": [[(-1, -1), (-1, 1), (1, -1), (1, 1)]],\n    \"E\": [[(1, 1), (-1, 1), (-1, -1), (1, -1)], [(-1, 0), (0.6, 0)]],\n    \"S\": [[(1, 1), (-1, 1), (-1, 0), (1, 0), (1, -1), (-1, -1)]],\n    \"W\": [[(-1, 1), (-0.5, -1), (0, 0.25), (0.5, -1), (1, 1)]],\n}\n\n\ndef add_letter(letter, cx, cy, s=0.05):\n    for stroke in LETTER_STROKES[letter]:\n        if boundary_data:\n            boundary_data.append(None)\n        for lx, ly in stroke:\n            boundary_data.append((round(cx + lx * s, 4), round(cy + ly * s, 4)))\n\n\nadd_letter(\"N\", 0.0, 1.16)\nadd_letter(\"E\", 1.16, 0.0)\nadd_letter(\"S\", 0.0, -1.16)\nadd_letter(\"W\", -1.16, 0.0)\n\n# Style - colors aligned to series add-order (poles reuse the feature color) --\ncustom_style = Style(\n    background=PAGE_BG,\n    plot_background=PAGE_BG,\n    foreground=INK,\n    foreground_strong=INK,\n    foreground_subtle=INK_SOFT,\n    colors=(BEDDING, BEDDING, FAULT, FAULT, JOINT, JOINT, DENSITY, INK_MUTED, INK),\n    title_font_size=66,\n    label_font_size=56,\n    major_label_font_size=44,\n    legend_font_size=44,\n    value_font_size=36,\n    opacity=0.78,\n    opacity_hover=0.95,\n)\n\nchart = pygal.XY(\n    width=2400,\n    height=2400,\n    style=custom_style,\n    title=\"stereonet-equal-area · python · pygal · anyplot.ai\",\n    show_legend=True,\n    legend_at_bottom=True,\n    legend_at_bottom_columns=4,\n    show_x_labels=False,\n    show_y_labels=False,\n    show_x_guides=False,\n    show_y_guides=False,\n    xrange=(-1.28, 1.28),\n    range=(-1.28, 1.28),\n    dots_size=0,\n    allow_interruptions=True,\n    margin_top=20,\n    margin_bottom=30,\n    margin_left=10,\n    margin_right=10,\n)\n\n# Great circles (planes) - legend entry per feature; line weight = hierarchy\n# Poles (dots) reuse the feature color but carry title=None so the legend stays\n# to four clean keys (Bedding, Faults, Joints, Pole density).\nfor name in (\"Bedding\", \"Faults\", \"Joints\"):\n    color = {\"Bedding\": BEDDING, \"Faults\": FAULT, \"Joints\": JOINT}[name]\n    chart.add(name, gc_series[name], stroke=True, show_dots=False, stroke_style={\"width\": gc_widths[name]})\n    chart.add(None, pole_series[name], stroke=False, show_dots=True, dots_size=pole_sizes[name])\n\n# Pole density contours - solid red, distinct from the dotted muted grid\nchart.add(\"Pole density\", contour_data, stroke=True, show_dots=False, stroke_style={\"width\": 2.2})\n\n# Equal-area net grid (hidden from legend) - light, dotted, recedes\nchart.add(None, grid_data, stroke=True, show_dots=False, stroke_style={\"width\": 1.0, \"dasharray\": \"1,6\"})\n\n# Primitive circle + ticks + cardinal letters (hidden from legend)\nchart.add(None, boundary_data, stroke=True, show_dots=False, stroke_style={\"width\": 2.6})\n\n# Save (theme-suffixed PNG + interactive HTML) -------------------------------\nchart.render_to_png(f\"plot-{THEME}.png\")\nwith open(f\"plot-{THEME}.html\", \"wb\") as f:\n    f.write(chart.render())\n"}