{"spec_id":"stereonet-equal-area","library":"bokeh","language":"python","code":"\"\"\" anyplot.ai\nstereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\nLibrary: bokeh 3.9.1 | Python 3.13.13\nQuality: 90/100 | Updated: 2026-06-16\n\"\"\"\n\nimport os\nimport sys\n\n\n# Prevent self-import: this file is named bokeh.py, which shadows the installed\n# bokeh package when its directory sits at the front of sys.path.\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 time\nfrom pathlib import Path\n\nimport numpy as np\nfrom bokeh.io import output_file, save\nfrom bokeh.models import ColorBar, ColumnDataSource, HoverTool, Label, Legend, LegendItem, LinearColorMapper\nfrom bokeh.plotting import figure\nfrom scipy.ndimage import gaussian_filter\nfrom scipy.stats import gaussian_kde\nfrom selenium import webdriver\nfrom selenium.webdriver.chrome.options import Options\n\n\n# Theme-adaptive chrome\nTHEME = os.getenv(\"ANYPLOT_THEME\", \"light\")\nPAGE_BG = \"#FAF8F1\" if THEME == \"light\" else \"#1A1A17\"\nELEVATED_BG = \"#FFFDF6\" if THEME == \"light\" else \"#242420\"\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 — Bedding/Joints/Faults are abstract geological\n# categories, so canonical order: green, lavender, blue.\nIMPRINT_PALETTE = [\"#009E73\", \"#C475FD\", \"#4467A3\", \"#BD8233\", \"#AE3030\", \"#2ABCCD\", \"#954477\", \"#99B314\"]\ncolors_map = {\"Bedding\": IMPRINT_PALETTE[0], \"Joints\": IMPRINT_PALETTE[1], \"Faults\": IMPRINT_PALETTE[2]}\n\n# Data - Synthetic structural geology measurements (strike, dip, feature_type)\nnp.random.seed(42)\n\n# Bedding planes: consistent NE strike with moderate dip\nbedding_strike = np.random.normal(45, 12, 40) % 360\nbedding_dip = np.random.normal(35, 8, 40).clip(5, 85)\n\n# Joint set: roughly E-W strike, steep dip\njoints_strike = np.random.normal(270, 15, 35) % 360\njoints_dip = np.random.normal(75, 10, 35).clip(5, 89)\n\n# Fault set: NW strike, moderate-steep dip\nfaults_strike = np.random.normal(315, 10, 25) % 360\nfaults_dip = np.random.normal(60, 12, 25).clip(5, 89)\n\nall_strikes = np.concatenate([bedding_strike, joints_strike, faults_strike])\nall_dips = np.concatenate([bedding_dip, joints_dip, faults_dip])\nall_types = [\"Bedding\"] * 40 + [\"Joints\"] * 35 + [\"Faults\"] * 25\n\n# Equal-area projection: convert pole (plunge, trend) to x, y\nR_net = 1.0\n\npole_trends = (all_strikes + 90 + 180) % 360\npole_plunges = 90.0 - all_dips\npole_trends_rad = np.radians(pole_trends)\npole_plunges_rad = np.radians(pole_plunges)\n\npole_r = R_net * np.sqrt(2) * np.sin((np.pi / 2 - pole_plunges_rad) / 2)\npole_x = pole_r * np.sin(pole_trends_rad)\npole_y = pole_r * np.cos(pole_trends_rad)\n\n# Great circles for each plane\ngc_xs = []\ngc_ys = []\ngc_types = []\n\nfor i in range(len(all_strikes)):\n    strike_rad = np.radians(all_strikes[i])\n    dip_rad = np.radians(all_dips[i])\n    dd_rad = strike_rad + np.pi / 2\n\n    sx = np.sin(strike_rad)\n    sy = np.cos(strike_rad)\n\n    dx = np.sin(dd_rad) * np.cos(dip_rad)\n    dy = np.cos(dd_rad) * np.cos(dip_rad)\n    dz = -np.sin(dip_rad)\n\n    alpha = np.linspace(0, np.pi, 90)\n    vx = np.cos(alpha) * sx + np.sin(alpha) * dx\n    vy = np.cos(alpha) * sy + np.sin(alpha) * dy\n    vz = np.sin(alpha) * dz\n\n    horiz = np.sqrt(vx**2 + vy**2)\n    plunge = np.arctan2(-vz, horiz)\n    trend = np.arctan2(vx, vy)\n\n    r = R_net * np.sqrt(2) * np.sin((np.pi / 2 - plunge) / 2)\n    gx = r * np.sin(trend)\n    gy = r * np.cos(trend)\n\n    # Clip great circle points near the primitive circle to reduce edge clutter\n    gc_dist = np.sqrt(gx**2 + gy**2)\n    keep = gc_dist <= 0.94 * R_net\n    gc_xs.append(gx[keep].tolist())\n    gc_ys.append(gy[keep].tolist())\n    gc_types.append(all_types[i])\n\n# Density grid for pole data using KDE\ngrid_n = 300\ngx_lin = np.linspace(-R_net, R_net, grid_n)\ngy_lin = np.linspace(-R_net, R_net, grid_n)\ngx_grid, gy_grid = np.meshgrid(gx_lin, gy_lin)\n\ndist_grid = np.sqrt(gx_grid**2 + gy_grid**2)\nmask = dist_grid <= R_net\n\npole_xy = np.vstack([pole_x, pole_y])\nkde = gaussian_kde(pole_xy, bw_method=0.2)\ndensity = kde(np.vstack([gx_grid.ravel(), gy_grid.ravel()])).reshape(grid_n, grid_n)\ndensity = gaussian_filter(density, sigma=3)\n\n# Normalize density inside the primitive circle\ndensity_masked = density.copy()\ndensity_masked[~mask] = np.nan\nd_min = np.nanmin(density_masked[mask])\nd_max = np.nanmax(density_masked[mask])\nd_norm = (density_masked - d_min) / (d_max - d_min)\n\n# Continuous data → Imprint sequential cmap (imprint_seq: brand green → blue).\nSEQ_LO = (0x00, 0x9E, 0x73)  # #009E73 brand green\nSEQ_HI = (0x44, 0x67, 0xA3)  # #4467A3 blue\nimprint_seq256 = [\n    \"#{:02X}{:02X}{:02X}\".format(\n        int(round(SEQ_LO[0] + (SEQ_HI[0] - SEQ_LO[0]) * t / 255.0)),\n        int(round(SEQ_LO[1] + (SEQ_HI[1] - SEQ_LO[1]) * t / 255.0)),\n        int(round(SEQ_LO[2] + (SEQ_HI[2] - SEQ_LO[2]) * t / 255.0)),\n    )\n    for t in range(256)\n]\n\n# Build uint32 RGBA density overlay from imprint_seq (vectorized)\nimg = np.zeros((grid_n, grid_n), dtype=np.uint32)\nview = img.view(dtype=np.uint8).reshape((grid_n, grid_n, 4))\nvisible = mask & (d_norm > 0.06)\nv = np.clip(np.where(visible, d_norm, 0.0), 0.0, 1.0)\nview[visible, 0] = (SEQ_LO[0] + (SEQ_HI[0] - SEQ_LO[0]) * v[visible]).astype(np.uint8)\nview[visible, 1] = (SEQ_LO[1] + (SEQ_HI[1] - SEQ_LO[1]) * v[visible]).astype(np.uint8)\nview[visible, 2] = (SEQ_LO[2] + (SEQ_HI[2] - SEQ_LO[2]) * v[visible]).astype(np.uint8)\nview[visible, 3] = (45 + 195 * v[visible]).astype(np.uint8)  # alpha rises with density\n\n# Plot - Square format for circular stereonet (2400 x 2400)\np = figure(\n    width=2400,\n    height=2400,\n    title=\"stereonet-equal-area · bokeh · anyplot.ai\",\n    x_range=(-1.36, 1.36),\n    y_range=(-1.40, 1.42),\n    tools=\"pan,wheel_zoom,reset,save\",\n    toolbar_location=None,\n    match_aspect=True,\n    min_border=20,\n)\n\n# Density heatmap as Bokeh image_rgba (distinctive Bokeh raster overlay)\np.image_rgba(image=[img], x=-R_net, y=-R_net, dw=2 * R_net, dh=2 * R_net, level=\"image\")\n\n# Density colorbar (imprint_seq) — inset into the empty bottom-left corner so it\n# sits with the circular plot rather than floating off to the side.\ncolor_mapper = LinearColorMapper(palette=imprint_seq256, low=0.0, high=1.0)\ncolor_bar = ColorBar(\n    color_mapper=color_mapper,\n    location=\"bottom_left\",\n    orientation=\"vertical\",\n    title=\"Pole density\",\n    title_text_font_size=\"26pt\",\n    title_text_font_style=\"italic\",\n    title_text_color=INK_SOFT,\n    major_label_text_font_size=\"24pt\",\n    major_label_text_color=INK_SOFT,\n    label_standoff=14,\n    width=42,\n    height=620,\n    padding=18,\n    margin=18,\n    background_fill_color=ELEVATED_BG,\n    background_fill_alpha=0.85,\n    border_line_color=INK_SOFT,\n    border_line_alpha=0.6,\n    major_tick_line_color=INK_SOFT,\n    major_tick_line_width=2,\n)\np.add_layout(color_bar)\n\n# Equal-area net grid lines (small circles at 10° dip intervals) — subtle\nfor dip_angle in range(10, 90, 10):\n    dip_rad = np.radians(dip_angle)\n    grid_r = R_net * np.sqrt(2) * np.sin(dip_rad / 2)\n    theta = np.linspace(0, 2 * np.pi, 180)\n    gx = grid_r * np.cos(theta)\n    gy = grid_r * np.sin(theta)\n    p.line(gx, gy, line_color=INK, line_width=1.5, line_alpha=0.14)\n\n# Great circle grid lines at every 30° azimuth — subtle\nfor az_deg in range(0, 180, 30):\n    az_rad = np.radians(az_deg)\n    alpha = np.linspace(0, np.pi, 90)\n    vx = np.cos(alpha) * np.sin(az_rad)\n    vy = np.cos(alpha) * np.cos(az_rad)\n    vz = -np.sin(alpha)\n    horiz = np.sqrt(vx**2 + vy**2)\n    plunge = np.arctan2(-vz, horiz)\n    trend = np.arctan2(vx, vy)\n    r = R_net * np.sqrt(2) * np.sin((np.pi / 2 - plunge) / 2)\n    grid_gx = r * np.sin(trend)\n    grid_gy = r * np.cos(trend)\n    p.line(grid_gx, grid_gy, line_color=INK, line_width=1.5, line_alpha=0.14)\n\n# Primitive circle (outer boundary)\ntheta_circle = np.linspace(0, 2 * np.pi, 360)\ncircle_x = R_net * np.cos(theta_circle)\ncircle_y = R_net * np.sin(theta_circle)\np.line(circle_x, circle_y, line_color=INK, line_width=4)\n\n# Tick marks every 10 degrees around perimeter\nfor deg in range(0, 360, 10):\n    rad = np.radians(deg)\n    tick_len = 0.06 if deg % 30 == 0 else 0.04\n    x_inner = (1.0 - tick_len) * R_net * np.sin(rad)\n    y_inner = (1.0 - tick_len) * R_net * np.cos(rad)\n    x_outer = 1.04 * R_net * np.sin(rad)\n    y_outer = 1.04 * R_net * np.cos(rad)\n    lw = 4 if deg % 90 == 0 else (3 if deg % 30 == 0 else 2)\n    p.line([x_inner, x_outer], [y_inner, y_outer], line_color=INK, line_width=lw)\n\n# Degree labels every 30 degrees (skip cardinals — labelled separately)\nfor deg in range(0, 360, 30):\n    if deg % 90 == 0:\n        continue\n    rad = np.radians(deg)\n    lx = 1.13 * R_net * np.sin(rad)\n    ly = 1.13 * R_net * np.cos(rad)\n    p.add_layout(\n        Label(\n            x=lx,\n            y=ly,\n            text=f\"{deg}°\",\n            text_font_size=\"22pt\",\n            text_align=\"center\",\n            text_baseline=\"middle\",\n            text_color=INK_SOFT,\n        )\n    )\n\n# Cardinal direction labels\nfor deg, label in [(0, \"N\"), (90, \"E\"), (180, \"S\"), (270, \"W\")]:\n    rad = np.radians(deg)\n    lx = 1.19 * R_net * np.sin(rad)\n    ly = 1.19 * R_net * np.cos(rad)\n    fs = \"40pt\" if label == \"N\" else \"32pt\"\n    p.add_layout(\n        Label(\n            x=lx,\n            y=ly,\n            text=label,\n            text_font_size=fs,\n            text_font_style=\"bold\",\n            text_align=\"center\",\n            text_baseline=\"middle\",\n            text_color=INK,\n        )\n    )\n\n# Great circles by feature type\nrenderers_gc = {}\nfor ftype in [\"Bedding\", \"Joints\", \"Faults\"]:\n    idxs = [j for j, t in enumerate(gc_types) if t == ftype]\n    fxs = [gc_xs[j] for j in idxs]\n    fys = [gc_ys[j] for j in idxs]\n    r = p.multi_line(fxs, fys, line_color=colors_map[ftype], line_width=1.6, line_alpha=0.22)\n    renderers_gc[ftype] = r\n\n# Poles by feature type with HoverTool\nrenderers_pole = {}\nfor ftype in [\"Bedding\", \"Joints\", \"Faults\"]:\n    idxs = [j for j, t in enumerate(all_types) if t == ftype]\n    px = pole_x[idxs]\n    py = pole_y[idxs]\n    strikes = all_strikes[idxs]\n    dips = all_dips[idxs]\n    source = ColumnDataSource(\n        data={\"x\": px, \"y\": py, \"strike\": np.round(strikes, 1), \"dip\": np.round(dips, 1), \"type\": [ftype] * len(idxs)}\n    )\n    r = p.scatter(\n        \"x\", \"y\", source=source, size=24, color=colors_map[ftype], line_color=PAGE_BG, line_width=2.5, alpha=0.95\n    )\n    renderers_pole[ftype] = r\n\n# HoverTool for pole data (Bokeh distinctive feature)\nhover = HoverTool(\n    renderers=list(renderers_pole.values()),\n    tooltips=[(\"Type\", \"@type\"), (\"Strike\", \"@strike°\"), (\"Dip\", \"@dip°\")],\n    point_policy=\"snap_to_data\",\n)\np.add_tools(hover)\n\n# Interactive legend (click to hide/show — Bokeh distinctive feature)\nlegend_items = []\nfor ftype in [\"Bedding\", \"Joints\", \"Faults\"]:\n    legend_items.append(LegendItem(label=ftype, renderers=[renderers_gc[ftype], renderers_pole[ftype]]))\n\nlegend = Legend(items=legend_items, location=\"top_right\")\nlegend.label_text_font_size = \"26pt\"\nlegend.label_text_color = INK_SOFT\nlegend.glyph_height = 40\nlegend.glyph_width = 40\nlegend.spacing = 16\nlegend.background_fill_color = ELEVATED_BG\nlegend.background_fill_alpha = 0.92\nlegend.border_line_color = INK_SOFT\nlegend.border_line_alpha = 0.6\nlegend.border_line_width = 2\nlegend.padding = 24\nlegend.margin = 20\nlegend.click_policy = \"hide\"\np.add_layout(legend)\n\n# Style\np.title.text_font_size = \"40pt\"\np.title.align = \"center\"\np.title.text_color = INK\np.xaxis.visible = False\np.yaxis.visible = False\np.xgrid.visible = False\np.ygrid.visible = False\np.outline_line_color = None\np.background_fill_color = PAGE_BG\np.border_fill_color = PAGE_BG\n\n# Subtitle annotation\np.add_layout(\n    Label(\n        x=0,\n        y=-1.31,\n        text=\"Lower-hemisphere equal-area (Schmidt) projection · Click legend to toggle\",\n        text_font_size=\"27pt\",\n        text_align=\"center\",\n        text_color=INK_MUTED,\n        text_font_style=\"italic\",\n    )\n)\n\n# Save — interactive HTML + headless-Chrome screenshot (export_png is unreliable here)\noutput_file(f\"plot-{THEME}.html\")\nsave(p)\n\nW, H = 2400, 2400\nopts = Options()\nfor arg in (\n    \"--headless=new\",\n    \"--no-sandbox\",\n    \"--disable-dev-shm-usage\",\n    \"--disable-gpu\",\n    f\"--window-size={W},{H}\",\n    \"--hide-scrollbars\",\n):\n    opts.add_argument(arg)\ndriver = webdriver.Chrome(options=opts)\ndriver.set_window_size(W, H)\ndriver.get(f\"file://{Path(f'plot-{THEME}.html').resolve()}\")\ntime.sleep(2)  # let bokeh's JS render the canvas\n\n# Headless Chrome's viewport is smaller than the requested window (window chrome),\n# which shrinks the screenshot below the figure's height. Compensate so the inner\n# viewport is exactly W x H, then re-render.\ninner_w = driver.execute_script(\"return window.innerWidth\")\ninner_h = driver.execute_script(\"return window.innerHeight\")\ndriver.set_window_size(W + (W - inner_w), H + (H - inner_h))\ntime.sleep(2)\ndriver.save_screenshot(f\"plot-{THEME}.png\")\ndriver.quit()\n"}