{"spec_id":"ternary-density","library":"pygal","language":"python","code":"\"\"\" anyplot.ai\nternary-density: Ternary Density Plot\nLibrary: pygal 3.1.0 | Python 3.13.13\nQuality: 81/100 | Updated: 2026-05-19\n\"\"\"\n\nimport math\nimport os\n\nimport cairosvg\nimport numpy as np\nfrom scipy.stats import gaussian_kde\n\n\n# Theme tokens\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\nH = math.sqrt(3) / 2\n\nnp.random.seed(42)\n\n# Data — sediment composition clusters (sand/silt/clay)\nn1 = 300\nsand1 = np.random.beta(5, 2, n1) * 60 + 35\nsilt1 = np.random.beta(2, 3, n1) * (100 - sand1) * 0.6\nclay1 = 100 - sand1 - silt1\n\nn2 = 250\nsilt2 = np.random.beta(5, 2, n2) * 50 + 40\nsand2 = np.random.beta(2, 4, n2) * (100 - silt2) * 0.5\nclay2 = 100 - sand2 - silt2\n\nn3 = 250\nclay3 = np.random.beta(4, 2, n3) * 45 + 40\nsand3 = np.random.beta(2, 5, n3) * (100 - clay3) * 0.4\nsilt3 = 100 - sand3 - clay3\n\nsand = np.clip(np.concatenate([sand1, sand2, sand3]), 0, 100)\nsilt = np.clip(np.concatenate([silt1, silt2, silt3]), 0, 100)\nclay = np.clip(np.concatenate([clay1, clay2, clay3]), 0, 100)\ntotal_comp = sand + silt + clay\nsand, silt, clay = sand / total_comp * 100, silt / total_comp * 100, clay / total_comp * 100\n\n# Ternary → Cartesian (Top=Clay, Bottom-left=Sand, Bottom-right=Silt)\nx_data = 0.5 * (2 * silt / 100 + clay / 100)\ny_data = H * clay / 100\n\n# Density grid\ngrid_res = 100\nxi = np.linspace(0, 1, grid_res)\nyi = np.linspace(0, H, grid_res)\nXi, Yi = np.meshgrid(xi, yi)\n\n# Triangle mask\nCi = Yi / H\nBi = Xi - Ci / 2\nmask = ((1 - Bi - Ci) >= -0.001) & (Bi >= -0.001) & (Ci >= -0.001)\n\n# KDE via scipy (Silverman bandwidth)\nkde = gaussian_kde(np.vstack([x_data, y_data]))\nZ = kde(np.vstack([Xi.ravel(), Yi.ravel()])).reshape(Xi.shape)\nZ = np.where(mask, Z, 0)\nZ_norm = (Z - Z.min()) / (Z.max() - Z.min() + 1e-10)\n\n# Viridis color table (11 waypoints, R/G/B ∈ [0,1])\nviridis_table = np.array(\n [\n [0.267004, 0.004874, 0.329415],\n [0.282327, 0.140926, 0.457517],\n [0.253935, 0.265254, 0.529983],\n [0.206756, 0.371758, 0.553117],\n [0.163625, 0.471133, 0.558148],\n [0.127568, 0.566949, 0.550556],\n [0.134692, 0.658636, 0.517649],\n [0.266941, 0.748751, 0.440573],\n [0.477504, 0.821444, 0.318195],\n [0.741388, 0.873449, 0.149561],\n [0.993248, 0.906157, 0.143936],\n ]\n)\n\n# Vectorized viridis color interpolation for all grid cells\nv_idx = np.minimum((Z_norm * 10).astype(int), 9)\nv_frac = Z_norm * 10 - v_idx.astype(float)\ncell_r = viridis_table[v_idx, 0] * (1 - v_frac) + viridis_table[v_idx + 1, 0] * v_frac\ncell_g = viridis_table[v_idx, 1] * (1 - v_frac) + viridis_table[v_idx + 1, 1] * v_frac\ncell_b = viridis_table[v_idx, 2] * (1 - v_frac) + viridis_table[v_idx + 1, 2] * v_frac\n\n# SVG canvas dimensions and coordinate mapping\nchart_width = 3600\nchart_height = 3600\nmargin = 200\nplot_size = chart_width - 2 * margin\nx_min, x_max = -0.12, 1.12\ny_min, y_max = -0.15, H + 0.15\nx_range = x_max - x_min\ny_range = y_max - y_min\n\n# Precompute pixel positions for density cell centers\ncs_x = xi[1] - xi[0]\ncs_y = yi[1] - yi[0]\npx_c = margin + (xi[:-1] + cs_x / 2 - x_min) / x_range * plot_size\npy_c = margin + (y_max - (yi[:-1] + cs_y / 2)) / y_range * plot_size\npw = cs_x / x_range * plot_size * 1.1\nph = cs_y / y_range * plot_size * 1.1\n\n# Density heatmap\ndensity_parts = ['\\n']\nfor i in range(grid_res - 1):\n for j in range(grid_res - 1):\n if mask[i, j] and Z_norm[i, j] > 0.05:\n r, g, b = int(cell_r[i, j] * 255), int(cell_g[i, j] * 255), int(cell_b[i, j] * 255)\n density_parts.append(\n f' \\n'\n )\ndensity_parts.append(\"\\n\")\ndensity_svg = \"\".join(density_parts)\n\n# Contour circles at density levels\ncontour_parts = [f'\\n']\nfor level in [0.2, 0.4, 0.6, 0.8]:\n thresh = Z_norm >= level\n for i in range(1, grid_res - 1):\n for j in range(1, grid_res - 1):\n if (\n thresh[i, j]\n and mask[i, j]\n and not all([thresh[i - 1, j], thresh[i + 1, j], thresh[i, j - 1], thresh[i, j + 1]])\n ):\n cx = margin + (xi[j] - x_min) / x_range * plot_size\n cy = margin + (y_max - yi[i]) / y_range * plot_size\n contour_parts.append(f' \\n')\ncontour_parts.append(\"\\n\")\ncontour_svg = \"\".join(contour_parts)\n\n# Triangle vertex pixel coordinates\nvx_sand = margin + (0.0 - x_min) / x_range * plot_size\nvy_sand = margin + (y_max - 0.0) / y_range * plot_size\nvx_silt = margin + (1.0 - x_min) / x_range * plot_size\nvy_silt = margin + (y_max - 0.0) / y_range * plot_size\nvx_clay = margin + (0.5 - x_min) / x_range * plot_size\nvy_clay = margin + (y_max - H) / y_range * plot_size\n\ntriangle_svg = f''\n\n# Grid lines at 25%, 50%, 75% (sparser spacing for cleaner readability)\ngrid_parts = [f'\\n']\nfor pct in [0.25, 0.50, 0.75]:\n # Constant clay lines: from (sand=1-pct, silt=0, clay=pct) to (sand=0, silt=1-pct, clay=pct)\n x1, y1 = 0.5 * pct, H * pct\n x2, y2 = 0.5 * (2 - pct), H * pct\n grid_parts.append(\n f' \\n'\n )\n # Constant silt lines: from (sand=1-pct, silt=pct, clay=0) to (sand=0, silt=pct, clay=1-pct)\n x1, y1 = pct, 0.0\n x2, y2 = 0.5 * (pct + 1), H * (1 - pct)\n grid_parts.append(\n f' \\n'\n )\n # Constant sand lines: from (sand=pct, silt=0, clay=1-pct) to (sand=pct, silt=1-pct, clay=0)\n x1, y1 = 0.5 * (1 - pct), H * (1 - pct)\n x2, y2 = 1 - pct, 0.0\n grid_parts.append(\n f' \\n'\n )\ngrid_parts.append(\"\\n\")\ngrid_svg = \"\".join(grid_parts)\n\n# Vertex labels with % unit indicator\nlabel_svg = f'''\n SAND (%)\n SILT (%)\n CLAY (%)\n\n'''\n\n# Percentage labels along edges at 25%, 50%, 75%\npct_parts = [f'\\n']\nfor pct in [25, 50, 75]:\n frac = pct / 100\n # Bottom edge (Sand-Silt axis): silt varies\n bx, by = frac, 0.0\n bpx = margin + (bx - x_min) / x_range * plot_size\n bpy = margin + (y_max - by) / y_range * plot_size\n pct_parts.append(f' {pct}%\\n')\n # Left edge (Sand-Clay axis): clay varies\n lx, ly = 0.5 * frac, H * frac\n lpx = margin + (lx - x_min) / x_range * plot_size\n lpy = margin + (y_max - ly) / y_range * plot_size\n pct_parts.append(f' {pct}%\\n')\n # Right edge (Silt-Clay axis): clay varies\n rx, ry = 0.5 * (2 - frac), H * frac\n rpx = margin + (rx - x_min) / x_range * plot_size\n rpy = margin + (y_max - ry) / y_range * plot_size\n pct_parts.append(f' {pct}%\\n')\npct_parts.append(\"\\n\")\npct_svg = \"\".join(pct_parts)\n\n# Colorbar\ncbar_x = chart_width - margin + 40\ncbar_y = margin + 200\ncbar_w = 40\ncbar_h = plot_size - 400\n\ncbar_parts = [\n '\\n',\n f' Density\\n',\n]\nn_steps = 50\nfor i in range(n_steps):\n val = 1 - i / n_steps\n vi = min(int(val * 10), 9)\n vf = val * 10 - vi\n cr = int((viridis_table[vi, 0] * (1 - vf) + viridis_table[vi + 1, 0] * vf) * 255)\n cg = int((viridis_table[vi, 1] * (1 - vf) + viridis_table[vi + 1, 1] * vf) * 255)\n cb = int((viridis_table[vi, 2] * (1 - vf) + viridis_table[vi + 1, 2] * vf) * 255)\n cy = cbar_y + i / n_steps * cbar_h\n ch = cbar_h / n_steps + 1\n cbar_parts.append(\n f' \\n'\n )\ncbar_parts += [\n f' \\n',\n f' High\\n',\n f' Low\\n',\n \"\\n\",\n]\ncolorbar_svg = \"\".join(cbar_parts)\n\n# Title\ntitle_svg = f'''\n Sediment Composition Distribution\n ternary-density · python · pygal · anyplot.ai\n\n'''\n\n# Clip path for density layer\nclip_svg = f'''\n \n \n \n'''\n\n# Assemble SVG\nsvg = f\"\"\"\n\"\"\"\n\n# Save HTML and PNG (theme-suffixed)\nwith open(f\"plot-{THEME}.html\", \"w\") as f:\n f.write(svg)\n\ncairosvg.svg2png(bytestring=svg.encode(\"utf-8\"), write_to=f\"plot-{THEME}.png\")\n"}