From 00551747579d4f19ba035072e49b219ff70c14b6 Mon Sep 17 00:00:00 2001 From: gronbow <59046437+gronbow@users.noreply.github.com> Date: Thu, 30 Jul 2026 09:29:37 +0800 Subject: [PATCH 1/2] Add Harker and TAS diagrams --- AGENTS.md | 4 +- README.md | 40 +- skills/geoskills/SKILL.md | 44 +- skills/geoskills/agents/openai.yaml | 4 +- .../classification/tas-lemaitre-2002.json | 156 ++++ .../examples/synthetic_major_element_data.csv | 11 + .../references/major-elements-and-tas.md | 52 ++ .../geoskills/scripts/inspect_major_data.py | 714 +++++++++++++++ .../geoskills/scripts/plot_geochem_common.py | 239 +++++ skills/geoskills/scripts/plot_harker.py | 770 ++++++++++++++++ skills/geoskills/scripts/plot_tas.py | 863 ++++++++++++++++++ tests/test_inspect_major_data.py | 146 +++ tests/test_plot_harker.py | 211 +++++ tests/test_plot_tas.py | 178 ++++ 14 files changed, 3417 insertions(+), 15 deletions(-) create mode 100644 skills/geoskills/assets/classification/tas-lemaitre-2002.json create mode 100644 skills/geoskills/examples/synthetic_major_element_data.csv create mode 100644 skills/geoskills/references/major-elements-and-tas.md create mode 100644 skills/geoskills/scripts/inspect_major_data.py create mode 100644 skills/geoskills/scripts/plot_geochem_common.py create mode 100644 skills/geoskills/scripts/plot_harker.py create mode 100644 skills/geoskills/scripts/plot_tas.py create mode 100644 tests/test_inspect_major_data.py create mode 100644 tests/test_plot_harker.py create mode 100644 tests/test_plot_tas.py diff --git a/AGENTS.md b/AGENTS.md index 99a9771..146977c 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -9,4 +9,6 @@ - Write machine-readable results to standard output and concise diagnostics to standard error where practical. - Add tests for numerical rules, invalid input, missing values, non-positive values, and exported files. - Do not claim the skill is research-ready until scientific fixtures and end-to-end tests pass. -- Preserve the published v0.1.0 REE workflow. GeoSkills v0.2.0 may add only trace-element spider diagrams; defer TAS, Harker, and other diagram families to later reviewed versions. +- Preserve the published v0.1.0 REE workflow and the reviewed v0.2.0 spider workflow. +- GeoSkills v0.3.0 may add Harker variation diagrams and volcanic TAS classification only; defer isotope, tectonic-discrimination, and other diagram families to later reviewed versions. +- Require explicit volcanic applicability and composition-basis declarations before TAS plotting; never silently classify intrusive or otherwise out-of-scope samples. diff --git a/README.md b/README.md index 0a377a9..1ea1570 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ GeoSkills 是一个面向地质学与地球化学研究的本地 Agent Skill。 公开的 v0.1.0 只提供球粒陨石标准化 REE 配分图。 -`feature/spider-diagram-v0.2` 分支是已完成人工审核的 v0.2.0 发布候选版：保留 REE，并新增原始地幔或 N-MORB 标准化微量元素蛛网图。TAS、Harker 和其他图解仍不在本版本范围内。 +v0.2.0 蛛网图工作流已经完成人工审核。当前本地 `feature/harker-tas-v0.3` 分支是 v0.3.0 审核候选版：保留前两类图件，并新增 Harker 变化图和带适用性保护的火山岩 TAS 分类图。该分支在用户审核前不会上传 GitHub。 ## 当前可以做什么 @@ -14,14 +14,16 @@ GeoSkills 是一个面向地质学与地球化学研究的本地 Agent Skill。 - 使用 Sun & McDonough（1989）C1 球粒陨石值标准化 La–Lu； - 使用 Sun & McDonough（1989）原始地幔、脚注明确修改的原始地幔或 N-MORB 值生成微量元素蛛网图； - 对明确标为 wt% 的 `K2O`、`P2O5`、`TiO2` 作可追溯的元素 ppm 换算； -- 生成投稿尺寸的 REE 配分图或蛛网图，并输出 SVG、PDF、600 dpi TIFF、600 dpi PNG； +- 以 SiO2 或用户指定变量为横轴，一次生成一幅或多幅 Harker 变化图； +- 使用 SiO2 与 Na2O + K2O 绘制火山岩 TAS 图，并输出逐样品分类与边界复核状态； +- 生成投稿尺寸的 REE、蛛网图、Harker 或 TAS 图件，并输出 SVG、PDF、600 dpi TIFF、600 dpi PNG； - 同步输出实际绘图数据 CSV 和机器可读的 JSON 运行报告； - 可选择完整四边框，以及经过碰撞检查的图内图例； - 全程在本地处理数据，绘图脚本不请求网络服务。 ## 安装到 Codex -以下公开仓库的默认分支目前仍是经过审核的 v0.1.0 REE 版本；v0.2.0 蛛网图发布候选版通过功能分支和 Pull Request 接受合并检查： +以下公开仓库的默认分支目前仍是经过审核的 v0.1.0 REE 版本；v0.2.0 蛛网图通过功能分支和 Pull Request 接受合并检查。v0.3.0 的 Harker/TAS 代码仍处于本地审核阶段： ```text https://github.com/gronbow/GeoSkills @@ -74,6 +76,29 @@ py -m venv .venv 默认方案为原文脚注推荐的 `pm-sm89-modified`；也可以明确选择 Table 1 印刷版 `pm-sm89` 或 `nmorb-sm89`。程序不会根据曲线形状替用户猜测标准化方案。 +检查 Harker/TAS 所需的主量元素数据： + +```powershell +.\.venv\Scripts\python.exe skills\geoskills\scripts\inspect_major_data.py skills\geoskills\examples\synthetic_major_element_data.csv +``` + +生成默认八面板 Harker 图： + +```powershell +.\.venv\Scripts\python.exe skills\geoskills\scripts\plot_harker.py skills\geoskills\examples\synthetic_major_element_data.csv --x SiO2 --y TiO2,Al2O3,Fe2O3T,MgO,CaO,Na2O,K2O,P2O5 --output-dir outputs\harker_figure +``` + +`--x` 和 `--y` 都可以按数据列自定义。程序默认不画回归线，也不会把相关性直接解释为分离结晶或岩浆混合。 +默认八面板会自动排成紧凑的 4 × 2 网格，并共用一个横轴标题；较长的分组图例会自动换行，避免超出图幅。 + +生成火山岩 TAS 图： + +```powershell +.\.venv\Scripts\python.exe skills\geoskills\scripts\plot_tas.py skills\geoskills\examples\synthetic_major_element_data.csv --confirm-volcanic --composition-basis anhydrous-normalized --output-dir outputs\tas_figure +``` + +TAS 命令要求明确确认样品属于火山岩，并声明数据是无水归一化值还是原始报告值。程序不会把侵入岩、碳酸岩或其他不适用样品静默套入火山岩名称。 + 程序不会修改原始表格。重复使用同一输出名称时，只有显式加入 `--overwrite` 才会替换已有结果。 ## 项目结构 @@ -87,7 +112,9 @@ GeoSkills/ │ └── geoskills/ │ ├── SKILL.md │ ├── agents/openai.yaml -│ ├── assets/normalization/ +│ ├── assets/ +│ │ ├── normalization/ +│ │ └── classification/ │ ├── examples/ │ ├── references/ │ └── scripts/ @@ -100,12 +127,15 @@ GeoSkills/ - 缺失值和低于检出限状态保留为曲线断点；零和负数会阻止对数坐标绘图。 - 标准化表包含文献、DOI、表格位置、版本和核对记录。 - 原始地幔表保留文献 Table 1 的 Cs、Pb 数值，并明确提示原文脚注中的 modified 版本；程序不会静默混用。 +- Harker 图只展示变量间的协变关系；默认不添加拟合线，也不从相关性单独推断岩浆过程。 +- TAS 边界使用带文献与版本信息的本地资产。恰好位于边界上的样品标为 `review_required`，不会静默选择一侧。 +- TAS 仅用于明确确认的火山岩；`as-reported` 结果标为初步分类，不能替代无水归一化后的专业判断。 - 图形可以展示富集程度、斜率和平行性，但不能单独证明岩浆源区、部分熔融、分离结晶或构造环境。 - `local_data/` 和 `outputs/` 已排除在 Git 之外；不要提交私人或未发表数据。 ## 当前状态 -GeoSkills v0.1.0 是仅包含 REE 配分图的公开测试版。v0.2.0 蛛网图发布候选版已通过全部 41 项自动测试、真实数据测试、导出审计和人工图形复核；远端 CI 与 Pull Request 合并完成后再作为稳定版本发布。 +GeoSkills v0.1.0 是仅包含 REE 配分图的公开测试版。v0.2.0 蛛网图发布候选版已通过 41 项自动测试、真实数据测试、导出审计和人工图形复核。v0.3.0 Harker/TAS 候选版已在本地通过 69 项自动测试、已发表数据的 Harker 验证、合成火山岩 TAS 验证、导出审计和人工图面复核；用户确认前不会提交或上传。 ## 许可 diff --git a/skills/geoskills/SKILL.md b/skills/geoskills/SKILL.md index 419ec55..cb08433 100644 --- a/skills/geoskills/SKILL.md +++ b/skills/geoskills/SKILL.md @@ -1,20 +1,22 @@ --- name: geoskills -description: "Inspect whole-rock geochemical tables and create validated, submission-oriented normalized element-pattern figures from CSV, TXT, or Excel data. Use GeoSkills for chondrite-normalized rare-earth-element patterns or primitive-mantle/N-MORB-normalized trace-element spider diagrams: identify sample, group, element, and supported oxide columns; verify units and invalid values; apply versioned Sun and McDonough (1989) references; export editable SVG/PDF plus 600 dpi TIFF/PNG and source data; or explain why a table cannot yet be plotted safely. Do not use the current version for TAS, Harker, isotope, discrimination, or other geochemical diagrams." +description: "Inspect whole-rock geochemical tables and create validated, submission-oriented REE patterns, trace-element spider diagrams, Harker variation diagrams, or volcanic TAS classification plots from CSV, TXT, or Excel data. Use GeoSkills to identify sample, group, element, and oxide columns; verify ppm or wt% units and invalid values; apply versioned Sun and McDonough (1989) references; require explicit TAS applicability declarations; export editable SVG/PDF plus 600 dpi TIFF/PNG and exact plotted data; or explain why a table cannot yet be plotted safely. Do not use the current version for isotope, tectonic-discrimination, or other diagram families." --- -# GeoSkills v0.2.0 release candidate +# GeoSkills v0.3.0 review candidate -Create reproducible REE patterns and trace-element spider diagrams through deterministic local Python scripts. Use the language model to select and explain the workflow, never to calculate normalized ratios manually. +Create reproducible REE patterns, trace-element spider diagrams, Harker variation diagrams, and volcanic TAS plots through deterministic local Python scripts. Use the language model to select and explain the workflow, never to calculate normalized ratios or classification fields manually. ## Development status -Treat v0.1.0 REE plotting as the stable published baseline. The v0.2.0 spider workflow is a reviewed release candidate: its scientific fixtures, real-data output, export QA, Skill validation, and user review have passed. Do not imply support for TAS, Harker, isotope, or discrimination diagrams. +Treat v0.1.0 REE plotting as the stable published baseline. The v0.2.0 spider workflow has passed scientific fixtures, real-data output, export QA, Skill validation, and user review. Harker and TAS are v0.3.0 review-candidate workflows and must not be described as published until their local review and release process are complete. Do not imply support for isotope or tectonic-discrimination diagrams. ## Route the request - For chondrite-normalized La–Lu patterns, use the REE workflow. - For multi-element primitive-mantle or N-MORB-normalized patterns, use the spider workflow. +- For one X variable against one or more oxides/elements, use the Harker workflow. +- For volcanic-rock classification by SiO2 and total alkalis, use the TAS workflow. - Stop and clarify when the requested diagram, unit, reference composition, sample column, or group column is ambiguous. ## REE workflow @@ -40,19 +42,42 @@ Use `--axes-frame full --legend-layout inside-auto` when a four-sided frame and The spider plot defaults to a full frame, no background grid, and `inside-auto` legend placement. It exports editable SVG/PDF, 600 dpi LZW TIFF, 600 dpi PNG, exact normalized source data, and a machine-readable report from the same figure object. +## Harker workflow + +1. Run `scripts/check_environment.py`. +2. Run `scripts/inspect_major_data.py INPUT`. +3. Confirm the worksheet, sample identifier, optional group, the X variable, requested Y variables, and each analyte's unit. +4. Run `scripts/plot_harker.py INPUT --x SiO2 --y TiO2,Al2O3,Fe2O3T,MgO,CaO,Na2O,K2O,P2O5 --output-dir OUTPUT_DIR`. +5. Inspect every panel at final size and return the SVG/PDF/TIFF/PNG, plotted-data CSV, and JSON report. + +`SiO2` is the conventional default X variable. The user may select another validated analyte explicitly. Do not add regression lines or assign a petrogenetic process from correlation alone. + +## TAS workflow + +1. Run `scripts/check_environment.py`. +2. Run `scripts/inspect_major_data.py INPUT`. +3. Confirm that the samples are volcanic and ask whether values are already anhydrous-normalized or should be plotted as reported. +4. Run `scripts/plot_tas.py INPUT --confirm-volcanic --composition-basis anhydrous-normalized --output-dir OUTPUT_DIR`. +5. Inspect the final-size figure, classification CSV, boundary/outside statuses, warnings, and JSON report. + +TAS uses `SiO2` versus `Na2O + K2O` in wt%. Never infer volcanic applicability from sample names. Do not use the volcanic fields for plutonic rocks, carbonatites, kimberlites, lamproites, or strongly altered compositions. Treat `as-reported` classifications as provisional and points on a boundary as `review_required`. + ## Scientific guardrails -- Require confirmed ppm units for direct element concentrations. +- Require explicit wt% units for major oxides and ppm units for direct elemental concentrations. +- Preserve blanks and below-detection-limit states as missing observations. Never replace them with zero or an invented detection limit. +- Reject negative concentrations and duplicate analyte mappings. Permit true zero only on linear Harker/TAS axes. - Convert only explicit `K2O`, `P2O5`, and `TiO2` wt% columns to K, P, and Ti ppm. Record the CIAAW/IUPAC atomic weights, formula, and factor used. -- Preserve blanks and below-detection-limit states as gaps. Never replace them with zero or an invented detection limit. - Reject finite zero and negative values on logarithmic axes. - Preserve the cited Sun and McDonough (1989) incompatibility order even when the user selects a subset. - Use only the versioned local assets `PrimitiveMantle_SM89`, `PrimitiveMantleModified_SM89`, and `NMORB_SM89` for the spider draft. - Use the source footnote's modified primitive mantle as the spider default. Keep the printed and modified variants separate; do not silently replace Cs or Pb, and warn when an affected element is plotted with the printed variant. - Set log limits from finite positive ratios, add a declared margin, round to clean decimal bounds without clipping, and show unity only when it lies inside the range. +- Use the versioned `TAS_LeMaitre2002_Volcanic_CombinedT` asset; keep `Trachyte/Trachydacite` and `Tephrite/Basanite` unresolved when normative information is absent. +- Round Harker limits outward to clean values without clipping data. Keep TAS at its declared fixed model limits so the field geometry is not distorted. - Use colour plus line style for groups and sample symbols so colour is not the sole identifier. - Keep user data local; plotting scripts must not make network requests. -- Describe enrichment, depletion, slopes, and visible anomalies conservatively. Do not assign a unique source, melting process, mineral control, alteration history, or tectonic setting from one normalized pattern. +- Describe enrichment, depletion, slopes, covariation, and visible anomalies conservatively. Do not assign a unique source, melting process, mineral control, alteration history, or tectonic setting from one plot. ## Related resources @@ -67,6 +92,11 @@ The spider plot defaults to a full frame, no background grid, and `inside-auto` | [assets/normalization/primitive-mantle-sm89.json](assets/normalization/primitive-mantle-sm89.json) | Audit primitive-mantle values and the Cs/Pb footnote | | [assets/normalization/primitive-mantle-modified-sm89.json](assets/normalization/primitive-mantle-modified-sm89.json) | Audit the explicit footnote-modified Cs/Pb variant | | [assets/normalization/nmorb-sm89.json](assets/normalization/nmorb-sm89.json) | Audit N-MORB values | +| [scripts/inspect_major_data.py](scripts/inspect_major_data.py) | Inspect major/trace analytes, units, missing states, and Harker/TAS readiness | +| [scripts/plot_harker.py](scripts/plot_harker.py) | Create customizable Harker grids | +| [scripts/plot_tas.py](scripts/plot_tas.py) | Create guarded volcanic TAS figures and classifications | +| [assets/classification/tas-lemaitre-2002.json](assets/classification/tas-lemaitre-2002.json) | Audit TAS fields, labels, applicability, and reference provenance | | [references/data-contract.md](references/data-contract.md) | Audit the REE input contract | | [references/scientific-method.md](references/scientific-method.md) | Audit the REE method | | [references/spider-method.md](references/spider-method.md) | Audit spider input, normalization, conversion, plotting, and interpretation rules | +| [references/major-elements-and-tas.md](references/major-elements-and-tas.md) | Audit Harker/TAS input, plotting, classification, and interpretation rules | diff --git a/skills/geoskills/agents/openai.yaml b/skills/geoskills/agents/openai.yaml index 29c7a29..1ad54af 100644 --- a/skills/geoskills/agents/openai.yaml +++ b/skills/geoskills/agents/openai.yaml @@ -1,4 +1,4 @@ interface: display_name: "GeoSkills" - short_description: "Create validated REE and trace-element pattern plots" - default_prompt: "Use $geoskills to inspect my geochemical table and create a validated REE or trace-element spider-plot figure bundle." + short_description: "Create validated REE, spider, Harker, and TAS plots" + default_prompt: "Use $geoskills to inspect my geochemical table and create a validated REE, trace-element spider, Harker, or volcanic TAS figure bundle." diff --git a/skills/geoskills/assets/classification/tas-lemaitre-2002.json b/skills/geoskills/assets/classification/tas-lemaitre-2002.json new file mode 100644 index 0000000..3060d5d --- /dev/null +++ b/skills/geoskills/assets/classification/tas-lemaitre-2002.json @@ -0,0 +1,156 @@ +{ + "schema_version": "1.0", + "id": "TAS_LeMaitre2002_Volcanic_CombinedT", + "display_name": "IUGS volcanic TAS classification", + "diagram_type": "total_alkali_silica", + "axes": { + "x": { + "analyte": "SiO2", + "unit": "wt%", + "limits": [35.0, 90.0] + }, + "y": { + "formula": "Na2O + K2O", + "unit": "wt%", + "limits": [0.0, 20.0] + } + }, + "applicability": { + "rock_domain": "volcanic", + "composition_basis": "volatile-free 100% basis recommended", + "not_for": [ + "plutonic-rock IUGS classification", + "carbonatites", + "kimberlites", + "lamproites", + "strongly altered compositions without independent review" + ] + }, + "source": { + "classification": "Le Maitre, R.W. (ed.) et al. (2002), Igneous Rocks: A Classification and Glossary of Terms, 2nd ed.", + "classification_doi": "10.1017/CBO9780511535581", + "boundary_construction": "Le Bas, M.J., Le Maitre, R.W. and Woolley, A.R. (1992), The construction of the Total Alkali-Silica chemical classification of volcanic rocks.", + "boundary_doi": "10.1007/BF01160698", + "note": "Volcanic names are used. The T field remains Trachyte/Trachydacite and U1 remains Tephrite/Basanite because those subtypes require additional normative information." + }, + "crosscheck": { + "project": "pyrolite", + "release": "v0.3.7", + "release_commit": "ab651b448ecd2ce1bc22d752b80a977a1f4dd1ca", + "path": "pyrolite/data/models/TAS/config_lemaitre_combined.json", + "verification_date": "2026-07-29", + "note": "Polygon vertices were checked against the pinned LeMaitreCombined configuration; volcanic field names and unresolved subtype labels are recorded explicitly here." + }, + "boundary_policy": { + "point_on_boundary": "review_required", + "tolerance_wt_percent": 1e-09, + "outside_model": "unclassified" + }, + "fields": [ + { + "id": "F", + "name": "Foidite", + "plot_label": "Foidite", + "label_position": [39.0, 11.0], + "polygon": [[35.0, 9.0], [37.0, 14.0], [52.5, 18.0], [52.5, 14.0], [48.4, 11.5], [45.0, 9.4], [41.0, 7.0], [41.0, 3.0], [37.0, 3.0]] + }, + { + "id": "Pc", + "name": "Picrobasalt", + "plot_label": "Picro-\nbasalt", + "label_position": [43.0, 1.4], + "polygon": [[41.0, 3.0], [45.0, 3.0], [45.0, 2.0], [45.0, 0.0], [41.0, 0.0]] + }, + { + "id": "B", + "name": "Basalt", + "plot_label": "Basalt", + "label_position": [48.5, 2.4], + "polygon": [[45.0, 0.0], [45.0, 5.0], [52.0, 5.0], [52.0, 0.0]] + }, + { + "id": "O1", + "name": "Basaltic andesite", + "plot_label": "Basaltic\nandesite", + "label_position": [54.5, 2.7], + "polygon": [[52.0, 0.0], [52.0, 5.0], [57.0, 5.9], [57.0, 0.0]] + }, + { + "id": "O2", + "name": "Andesite", + "plot_label": "Andesite", + "label_position": [60.0, 3.0], + "polygon": [[57.0, 0.0], [57.0, 5.9], [63.0, 7.0], [63.0, 0.0]] + }, + { + "id": "O3", + "name": "Dacite", + "plot_label": "Dacite", + "label_position": [66.2, 4.2], + "polygon": [[63.0, 0.0], [63.0, 7.0], [69.0, 8.0], [77.3, 0.0]] + }, + { + "id": "R", + "name": "Rhyolite", + "plot_label": "Rhyolite", + "label_position": [75.0, 9.3], + "polygon": [[69.0, 8.0], [69.0, 13.0], [85.9, 6.8], [87.5, 4.7], [77.3, 0.0]] + }, + { + "id": "S1", + "name": "Trachybasalt", + "plot_label": "Trachy-\nbasalt", + "label_position": [49.2, 5.8], + "polygon": [[45.0, 5.0], [49.4, 7.3], [52.0, 5.0]] + }, + { + "id": "S2", + "name": "Basaltic trachyandesite", + "plot_label": "Basaltic\ntrachyandesite", + "label_position": [53.5, 7.2], + "polygon": [[49.4, 7.3], [53.0, 9.3], [57.0, 5.9], [52.0, 5.0]] + }, + { + "id": "S3", + "name": "Trachyandesite", + "plot_label": "Trachy-\nandesite", + "label_position": [58.4, 8.5], + "polygon": [[53.0, 9.3], [57.6, 11.7], [61.0, 8.6], [63.0, 7.0], [57.0, 5.9]] + }, + { + "id": "T", + "name": "Trachyte/Trachydacite", + "plot_label": "Trachyte/\nTrachydacite", + "label_position": [65.2, 10.7], + "polygon": [[57.6, 11.7], [61.0, 13.5], [63.0, 16.2], [69.0, 13.0], [69.0, 8.0], [63.0, 7.0], [61.0, 8.6]] + }, + { + "id": "U1", + "name": "Tephrite/Basanite", + "plot_label": "Tephrite/\nBasanite", + "label_position": [44.0, 5.7], + "polygon": [[41.0, 3.0], [41.0, 7.0], [45.0, 9.4], [49.4, 7.3], [45.0, 5.0], [45.0, 3.0]] + }, + { + "id": "U2", + "name": "Phonotephrite", + "plot_label": "Phono-\ntephrite", + "label_position": [48.6, 9.4], + "polygon": [[45.0, 9.4], [48.4, 11.5], [53.0, 9.3], [49.4, 7.3]] + }, + { + "id": "U3", + "name": "Tephriphonolite", + "plot_label": "Tephri-\nphonolite", + "label_position": [53.4, 11.7], + "polygon": [[48.4, 11.5], [52.5, 14.0], [57.6, 11.7], [53.0, 9.3]] + }, + { + "id": "Ph", + "name": "Phonolite", + "plot_label": "Phonolite", + "label_position": [57.6, 14.5], + "polygon": [[52.5, 14.0], [52.5, 18.0], [57.0, 18.0], [63.0, 16.2], [61.0, 13.5], [57.6, 11.7]] + } + ] +} diff --git a/skills/geoskills/examples/synthetic_major_element_data.csv b/skills/geoskills/examples/synthetic_major_element_data.csv new file mode 100644 index 0000000..24c8a1c --- /dev/null +++ b/skills/geoskills/examples/synthetic_major_element_data.csv @@ -0,0 +1,11 @@ +Sample,Group,SiO2_wt%,TiO2_wt%,Al2O3_wt%,Fe2O3T_wt%,MnO_wt%,MgO_wt%,CaO_wt%,Na2O_wt%,K2O_wt%,P2O5_wt%,Rb_ppm,Sr_ppm,Zr_ppm +SYN-01,Series A,47.0,1.45,16.2,10.8,0.18,7.8,10.2,2.5,0.5,0.28,22,430,105 +SYN-02,Series A,50.0,1.20,16.8,9.4,0.16,5.9,8.6,3.5,2.5,0.24,48,390,145 +SYN-03,Series A,54.0,0.98,17.1,8.1,0.14,4.1,7.1,3.0,1.0,0.20,62,350,170 +SYN-04,Series B,54.0,0.92,17.5,7.9,0.13,3.7,6.8,4.2,3.0,0.18,85,320,205 +SYN-05,Series B,59.5,0.70,16.9,6.4,0.11,2.4,5.1,3.5,1.3,0.14,105,270,240 +SYN-06,Series B,59.0,0.68,17.0,6.2,0.10,2.2,4.8,5.0,4.0,0.13,132,245,275 +SYN-07,Series C,66.0,0.42,15.8,4.6,0.08,1.2,3.2,4.0,1.5,0.09,155,190,315 +SYN-08,Series C,74.0,0.20,13.9,2.5,0.05,0.4,1.4,4.0,2.2,0.04,190,115,365 +SYN-09,Series C,49.0,1.10,17.8,8.8,0.15,4.8,7.4,6.0,4.3,0.22,118,360,230 +SYN-10,Series C,58.0,0.55,18.1,5.5,0.09,1.8,3.8,7.0,7.0,0.11,175,205,340 diff --git a/skills/geoskills/references/major-elements-and-tas.md b/skills/geoskills/references/major-elements-and-tas.md new file mode 100644 index 0000000..ea7dc6a --- /dev/null +++ b/skills/geoskills/references/major-elements-and-tas.md @@ -0,0 +1,52 @@ +# Major-element, Harker, and TAS contract + +## Harker purpose + +Use Harker variation diagrams to show how selected oxides or elements covary with a declared differentiation index. Use `SiO2` as the conventional default x variable, but permit another validated analyte when the user explicitly selects it. + +Treat the panels as comparative evidence. Do not infer fractional crystallization, magma mixing, assimilation, source variation, or alteration from correlation alone. Do not add regression lines by default. + +## Accepted input + +- One unambiguous sample identifier and an optional group column. +- At least two explicitly unit-labelled geochemical analytes. +- Major oxides in wt% and direct elemental concentrations in ppm. +- Flat row-per-sample tables or unambiguous published-supplement tables with samples in columns and unit-labelled major/trace sections. + +Preserve blanks and below-detection-limit states as missing points. Reject non-numeric values, duplicate analyte mappings, and negative concentrations. Zero is permitted on linear axes. + +## TAS applicability + +TAS uses: + +```text +x = SiO2 (wt%) +y = Na2O (wt%) + K2O (wt%) +``` + +Require an explicit confirmation that the samples are volcanic. Require the user to state whether the input is already on a volatile-free/anhydrous 100% basis or is being plotted as reported. The current version records that declaration but does not silently recalculate analyses. + +Do not use volcanic TAS names as the IUGS classification of plutonic rocks, carbonatites, kimberlites, lamproites, or strongly altered compositions. Report an as-reported classification as provisional. + +## TAS fields + +Use the versioned `TAS_LeMaitre2002_Volcanic_CombinedT` asset. It records the volcanic field names recommended in Le Maitre et al. (2002), the boundary-construction paper by Le Bas et al. (1992), and a pinned pyrolite v0.3.7 cross-check. + +Keep `Trachyte/Trachydacite` and `Tephrite/Basanite` unresolved because those distinctions need normative information not supplied by SiO2 and total alkalis alone. Mark an analysis exactly on a field boundary as `review_required`; do not choose a side silently. + +## Plot and export + +- Use a white background and full four-sided axes by default. +- Use colour plus marker shape for groups. +- Use one shared, width-checked legend and one common X-axis title for a Harker grid. +- Choose a compact Harker panel layout automatically while retaining an explicit column override. +- Place the TAS legend inside only when it does not overlap points or field labels; otherwise move it outside-right. +- Use clean outward-rounded linear limits for Harker panels without clipping data. +- Keep the TAS model at its fixed declared limits so field geometry is not visually distorted. +- Export SVG, PDF, 600 dpi LZW TIFF, 600 dpi PNG, exact source-data CSV, and a JSON report from the same Matplotlib figure. + +## References + +- Harker, A. (1909), *The Natural History of Igneous Rocks*. +- Le Bas, M.J., Le Maitre, R.W. and Woolley, A.R. (1992), *Mineralogy and Petrology* 46, 1–22. DOI: `10.1007/BF01160698`. +- Le Maitre, R.W. (ed.) et al. (2002), *Igneous Rocks: A Classification and Glossary of Terms*, 2nd ed. DOI: `10.1017/CBO9780511535581`. diff --git a/skills/geoskills/scripts/inspect_major_data.py b/skills/geoskills/scripts/inspect_major_data.py new file mode 100644 index 0000000..7575ff5 --- /dev/null +++ b/skills/geoskills/scripts/inspect_major_data.py @@ -0,0 +1,714 @@ +#!/usr/bin/env python3 +"""Inspect major-element and Harker/TAS input without modifying the source.""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import re +import sys +import unicodedata +from pathlib import Path +from typing import Any + +import pandas as pd + +from inspect_data import ( + BDL_PATTERN, + GROUP_NAMES, + MAX_FILE_SIZE_BYTES, + SAMPLE_NAMES, + SUPPORTED_SUFFIXES, + InspectionError, + clean_name, + issue, + problem_examples, + resolve_sheet, + sniff_text_format, +) +from inspect_spider_data import SPIDER_ELEMENT_ORDER + + +MAJOR_OXIDE_ALIASES = { + "sio2": "SiO2", + "silica": "SiO2", + "tio2": "TiO2", + "al2o3": "Al2O3", + "fe2o3t": "Fe2O3T", + "tfe2o3": "Fe2O3T", + "fe2o3total": "Fe2O3T", + "fe2o3tot": "Fe2O3T", + "feot": "FeOT", + "tfeo": "FeOT", + "feototal": "FeOT", + "feotot": "FeOT", + "fe2o3": "Fe2O3", + "feo": "FeO", + "mno": "MnO", + "mgo": "MgO", + "cao": "CaO", + "na2o": "Na2O", + "k2o": "K2O", + "p2o5": "P2O5", + "cr2o3": "Cr2O3", + "nio": "NiO", + "h2o": "H2O", + "h2oplus": "H2O+", + "h2ominus": "H2O-", + "co2": "CO2", + "loi": "LOI", + "lossonignition": "LOI", + "total": "Total", +} +MAJOR_OXIDE_ORDER = [ + "SiO2", + "TiO2", + "Al2O3", + "Fe2O3T", + "FeOT", + "Fe2O3", + "FeO", + "MnO", + "MgO", + "CaO", + "Na2O", + "K2O", + "P2O5", + "Cr2O3", + "NiO", + "H2O", + "H2O+", + "H2O-", + "CO2", + "LOI", + "Total", +] +EXTRA_TRACE_ELEMENTS = [ + "Sc", + "V", + "Cr", + "Co", + "Ni", + "Cu", + "Zn", + "Ga", + "Ge", + "As", + "Se", + "Br", + "Ag", + "Cd", + "In", + "Te", +] +TRACE_ELEMENT_ORDER = list( + dict.fromkeys([*EXTRA_TRACE_ELEMENTS, *SPIDER_ELEMENT_ORDER]) +) +TRACE_BY_KEY = {clean_name(element): element for element in TRACE_ELEMENT_ORDER} +UNIT_SUFFIX_PATTERN = re.compile( + r"[\s_\-\(\[]*(?:ppm|ppb|wt\s*\.?\s*(?:%|pct|percent))" + r"[\s\)\]]*\s*$", + re.IGNORECASE, +) + + +def infer_geochem_unit(value: object) -> str: + """Infer a unit only when it is explicitly present.""" + normalized = unicodedata.normalize("NFKC", str(value)).lower() + if re.search(r"(?:^|[^a-z])ppm(?:$|[^a-z])", normalized): + return "ppm" + if re.search(r"(?:^|[^a-z])ppb(?:$|[^a-z])", normalized): + return "ppb" + if re.search(r"(?:^|[^a-z])wt\s*\.?\s*(?:%|pct|percent)", normalized): + return "wt%" + return "unknown" + + +def match_geochem_analyte(value: object) -> dict[str, str] | None: + """Map an exact oxide or common trace-element label to a canonical name.""" + normalized = unicodedata.normalize("NFKC", str(value)).strip() + base = UNIT_SUFFIX_PATTERN.sub("", normalized) + key = clean_name(base) + if key in MAJOR_OXIDE_ALIASES: + analyte = MAJOR_OXIDE_ALIASES[key] + return { + "analyte": analyte, + "kind": "major_oxide", + "required_unit": "wt%", + } + if key in TRACE_BY_KEY: + analyte = TRACE_BY_KEY[key] + return { + "analyte": analyte, + "kind": "trace_element", + "required_unit": "ppm", + } + return None + + +def adapt_transposed_major_table( + raw: pd.DataFrame, +) -> tuple[pd.DataFrame, dict[str, Any]] | None: + """Convert a common analytes-by-row supplement into one row per sample.""" + if raw.empty or raw.shape[0] < 7 or raw.shape[1] < 3: + return None + + best_label_column: int | None = None + best_rows: list[tuple[int, dict[str, str]]] = [] + for column_index in range(raw.shape[1]): + matches: list[tuple[int, dict[str, str]]] = [] + seen: set[str] = set() + duplicate = False + for row_index, value in raw.iloc[:, column_index].items(): + analyte = match_geochem_analyte(value) + if analyte is None: + continue + if analyte["analyte"] in seen: + duplicate = True + break + seen.add(analyte["analyte"]) + matches.append((int(row_index), analyte)) + if not duplicate and len(matches) > len(best_rows): + best_label_column = column_index + best_rows = matches + + if best_label_column is None or len(best_rows) < 3: + return None + + sample_rows = [ + int(row_index) + for row_index, value in raw.iloc[:, best_label_column].items() + if clean_name(value) in SAMPLE_NAMES + and int(raw.iloc[int(row_index)].notna().sum()) >= 3 + ] + if len(sample_rows) != 1: + return None + sample_row = sample_rows[0] + + sample_columns = [ + column_index + for column_index in range(raw.shape[1]) + if column_index != best_label_column + and pd.notna(raw.iat[sample_row, column_index]) + and str(raw.iat[sample_row, column_index]).strip() + ] + if len(sample_columns) < 2: + return None + + group_row: int | None = None + for row_index in range(sample_row): + if clean_name(raw.iat[row_index, best_label_column]) in GROUP_NAMES: + group_row = row_index + + unit_by_row: dict[int, str] = {} + current_unit = "unknown" + matched_rows = {row_index for row_index, _ in best_rows} + for row_index in range(sample_row + 1, raw.shape[0]): + label = raw.iat[row_index, best_label_column] + explicit_unit = infer_geochem_unit(label) + if explicit_unit != "unknown" and row_index not in matched_rows: + current_unit = explicit_unit + if row_index in matched_rows: + unit_by_row[row_index] = ( + explicit_unit if explicit_unit != "unknown" else current_unit + ) + + converted: dict[str, list[Any]] = { + "Sample": [ + str(raw.iat[sample_row, column]).strip() + for column in sample_columns + ] + } + if group_row is not None: + group_values = pd.Series( + [raw.iat[group_row, column] for column in sample_columns], + dtype="object", + ).ffill() + if group_values.notna().any(): + converted["Group"] = [ + None if pd.isna(value) else str(value).strip() + for value in group_values.tolist() + ] + + recognized: list[dict[str, Any]] = [] + output_names = set(converted) + for row_index, analyte in best_rows: + unit = unit_by_row.get(row_index, "unknown") + output_name = ( + f"{analyte['analyte']}_{unit}" + if unit != "unknown" + else analyte["analyte"] + ) + if output_name in output_names: + return None + output_names.add(output_name) + converted[output_name] = [ + raw.iat[row_index, column] for column in sample_columns + ] + recognized.append( + { + **analyte, + "source_row": row_index + 1, + "source_label": str( + raw.iat[row_index, best_label_column] + ), + "inferred_unit": unit, + "output_column": output_name, + } + ) + + return pd.DataFrame(converted), { + "method": "auto_transpose_geochemical_analytes_by_row", + "sample_header_row": sample_row + 1, + "sample_count": len(sample_columns), + "element_label_column": best_label_column + 1, + "group_header_row": None if group_row is None else group_row + 1, + "recognized_analytes": recognized, + } + + +def read_major_table( + path: Path, + requested_sheet: str | None, +) -> tuple[pd.DataFrame | None, dict[str, Any]]: + """Read one supported table with major-element transposed-layout handling.""" + if not path.exists() or not path.is_file(): + raise InspectionError("找不到输入文件。") + if path.stat().st_size > MAX_FILE_SIZE_BYTES: + raise InspectionError("文件超过 20 MB；当前版本暂不处理更大的文件。") + suffix = path.suffix.lower() + if suffix not in SUPPORTED_SUFFIXES: + raise InspectionError("不支持该文件格式；请使用 .csv、.txt 或 .xlsx。") + + source: dict[str, Any] = { + "file_sha256": hashlib.sha256(path.read_bytes()).hexdigest(), + "size_bytes": path.stat().st_size, + "format": suffix, + "sheet": None, + "sheet_names": [], + "encoding": None, + "delimiter": None, + "layout": None, + "transformation": None, + } + try: + if suffix in {".csv", ".txt"}: + encoding, delimiter = sniff_text_format(path) + source["encoding"] = encoding + source["delimiter"] = "TAB" if delimiter == "\t" else delimiter + raw = pd.read_csv( + path, sep=delimiter, encoding=encoding, header=None + ) + adapted = adapt_transposed_major_table(raw) + if adapted is not None: + frame, transformation = adapted + source["layout"] = "column_per_sample_transposed" + source["transformation"] = transformation + return frame, source + source["layout"] = "row_per_sample" + return pd.read_csv( + path, sep=delimiter, encoding=encoding + ), source + + workbook = pd.ExcelFile(path) + source["sheet_names"] = workbook.sheet_names + selected_sheet = resolve_sheet( + workbook.sheet_names, requested_sheet + ) + if selected_sheet is None: + return None, source + source["sheet"] = selected_sheet + raw = pd.read_excel( + workbook, sheet_name=selected_sheet, header=None + ) + adapted = adapt_transposed_major_table(raw) + if adapted is not None: + frame, transformation = adapted + source["layout"] = "column_per_sample_transposed" + source["transformation"] = transformation + return frame, source + source["layout"] = "row_per_sample" + return pd.read_excel( + workbook, sheet_name=selected_sheet + ), source + except InspectionError: + raise + except Exception as exc: + raise InspectionError(f"无法读取表格：{exc}") from exc + + +def inspect_major_frame( + frame: pd.DataFrame, + source: dict[str, Any], + requested_sample_column: str | None = None, + requested_group_column: str | None = None, +) -> dict[str, Any]: + """Inspect identifiers, analyte mappings, units, and cell states.""" + columns = [str(column) for column in frame.columns] + automatic_sample_candidates = [ + column for column in frame.columns + if clean_name(column) in SAMPLE_NAMES + ] + automatic_group_candidates = [ + column for column in frame.columns + if clean_name(column) in GROUP_NAMES + ] + requested_sample_matches = [ + column for column in frame.columns + if str(column) == requested_sample_column + ] + requested_group_matches = [ + column for column in frame.columns + if str(column) == requested_group_column + ] + sample_candidates = ( + requested_sample_matches + if requested_sample_column is not None + else automatic_sample_candidates + ) + group_candidates = ( + requested_group_matches + if requested_group_column is not None + else automatic_group_candidates + ) + sample_column = ( + sample_candidates[0] if len(sample_candidates) == 1 else None + ) + issues: list[dict[str, Any]] = [] + + if frame.empty: + issues.append(issue("E601", "error", "表格没有数据行。")) + if ( + requested_sample_column is not None + and len(requested_sample_matches) != 1 + ): + issues.append( + issue( + "E606", + "error", + "指定的样品编号列不存在或不唯一。", + requested=requested_sample_column, + ) + ) + elif not sample_candidates: + issues.append( + issue("E602", "review", "未自动识别样品编号列，请明确指定。") + ) + elif len(sample_candidates) > 1: + issues.append( + issue( + "E602", + "review", + "识别到多个可能的样品编号列，请确认使用哪一列。", + columns=[str(column) for column in sample_candidates], + ) + ) + if ( + requested_group_column is not None + and len(requested_group_matches) != 1 + ): + issues.append( + issue( + "E607", + "error", + "指定的分组列不存在或不唯一。", + requested=requested_group_column, + ) + ) + + if sample_column is not None: + sample_series = frame[sample_column] + sample_text = sample_series.astype("string").str.strip() + missing_mask = sample_series.isna() | sample_text.eq("") + duplicate_mask = ( + ~missing_mask & sample_text.duplicated(keep=False) + ) + if int(missing_mask.sum()): + issues.append( + issue( + "E604", + "review", + "样品编号列包含空白值。", + column=str(sample_column), + ) + ) + if int(duplicate_mask.sum()): + issues.append( + issue( + "E605", + "review", + "样品编号存在重复值，请确认这些行是否代表同一样品。", + column=str(sample_column), + values=sorted( + sample_text[duplicate_mask] + .dropna() + .unique() + .tolist() + )[:10], + ) + ) + + recognized: list[dict[str, Any]] = [] + by_analyte: dict[str, list[str]] = {} + for column in frame.columns: + analyte = match_geochem_analyte(column) + if analyte is None: + continue + by_analyte.setdefault(analyte["analyte"], []).append(str(column)) + series = frame[column] + text_values = series.astype("string").str.strip() + missing_mask = series.isna() | text_values.eq("") + bdl_mask = ( + text_values.str.match(BDL_PATTERN, na=False) & ~missing_mask + ) + numeric = pd.to_numeric( + series.where(~missing_mask & ~bdl_mask), + errors="coerce", + ) + nonnumeric_mask = ( + ~missing_mask & ~bdl_mask & numeric.isna() + ) + negative_mask = numeric.notna() & numeric.lt(0) + unit = infer_geochem_unit(column) + recognized.append( + { + **analyte, + "column": str(column), + "unit": unit, + "missing": int(missing_mask.sum()), + "below_detection_limit": int(bdl_mask.sum()), + "non_numeric": int(nonnumeric_mask.sum()), + "negative": int(negative_mask.sum()), + } + ) + + if unit != analyte["required_unit"]: + message = ( + f"{analyte['analyte']} 列名或所在分区没有明确单位" + f" {analyte['required_unit']}。" + if unit == "unknown" + else ( + f"{analyte['analyte']} 必须使用" + f" {analyte['required_unit']}，但当前标记为 {unit}。" + ) + ) + issues.append( + issue( + "E621", + "review", + message, + analyte=analyte["analyte"], + column=str(column), + ) + ) + if int(bdl_mask.sum()): + issues.append( + issue( + "W631", + "warning", + "检测到低于检出限的值；程序保留为缺失点，没有替换成零。", + analyte=analyte["analyte"], + column=str(column), + examples=problem_examples( + frame, bdl_mask, column, sample_column + ), + ) + ) + if int(nonnumeric_mask.sum()): + issues.append( + issue( + "E631", + "review", + "检测到无法转换为数字的单元格。", + analyte=analyte["analyte"], + column=str(column), + examples=problem_examples( + frame, nonnumeric_mask, column, sample_column + ), + ) + ) + if int(negative_mask.sum()): + issues.append( + issue( + "E632", + "review", + "检测到负浓度，不能安全用于主量元素或 Harker/TAS 图解。", + analyte=analyte["analyte"], + column=str(column), + examples=problem_examples( + frame, negative_mask, column, sample_column + ), + ) + ) + if int(missing_mask.sum()): + issues.append( + issue( + "W611", + "warning", + "检测到缺失值；绘图时跳过对应的散点。", + analyte=analyte["analyte"], + column=str(column), + count=int(missing_mask.sum()), + ) + ) + + if len({item["analyte"] for item in recognized}) < 2: + issues.append( + issue( + "E603", + "error", + "至少需要识别到两个明确标注单位的地球化学变量。", + ) + ) + duplicates = { + analyte: names + for analyte, names in by_analyte.items() + if len(names) > 1 + } + if duplicates: + issues.append( + issue( + "E611", + "review", + "同一变量对应多个输入列，请明确选择。", + columns=duplicates, + ) + ) + + severities = {item["severity"] for item in issues} + if "error" in severities: + status = "error" + elif "review" in severities: + status = "needs_review" + else: + status = "ready" + return { + "status": status, + "source": source, + "shape": { + "rows": int(frame.shape[0]), + "columns": int(frame.shape[1]), + }, + "columns": columns, + "sample_id_candidates": [ + str(column) for column in sample_candidates + ], + "group_candidates": [ + str(column) for column in group_candidates + ], + "analytes": { + "major_oxide_order": MAJOR_OXIDE_ORDER, + "trace_element_order": TRACE_ELEMENT_ORDER, + "recognized": recognized, + }, + "issues": issues, + } + + +def prepare_geochem_frame( + frame: pd.DataFrame, + inspection: dict[str, Any], +) -> pd.DataFrame: + """Create canonical numeric columns from a ready inspection.""" + if inspection["status"] != "ready": + raise InspectionError( + f"输入检查状态为 {inspection['status']}，不能安全准备绘图数据。" + ) + sample_candidates = inspection["sample_id_candidates"] + if len(sample_candidates) != 1: + raise InspectionError("必须明确识别一个样品编号列。") + output_columns = [sample_candidates[0]] + if len(inspection["group_candidates"]) == 1: + output_columns.append(inspection["group_candidates"][0]) + output = frame.loc[:, output_columns].copy() + for item in inspection["analytes"]["recognized"]: + output[item["analyte"]] = pd.to_numeric( + frame[item["column"]], errors="coerce" + ) + return output + + +def major_error_report(path: Path, message: str) -> dict[str, Any]: + return { + "status": "error", + "operation": "major_element_inspection", + "source": {"format": path.suffix.lower()}, + "issues": [issue("E600", "error", message)], + } + + +def inspect_major_path( + path: Path, + requested_sheet: str | None = None, + requested_sample_column: str | None = None, + requested_group_column: str | None = None, +) -> dict[str, Any]: + """Inspect one major-element input path and return a JSON-ready report.""" + try: + frame, source = read_major_table(path, requested_sheet) + if frame is None: + return { + "status": "needs_sheet", + "operation": "major_element_inspection", + "source": source, + "issues": [ + issue( + "E111", + "review", + "Excel 文件包含多个工作表，请使用 --sheet 指定一个工作表。", + sheet_names=source["sheet_names"], + ) + ], + } + report = inspect_major_frame( + frame, + source, + requested_sample_column, + requested_group_column, + ) + report["operation"] = "major_element_inspection" + return report + except (InspectionError, OSError) as exc: + return major_error_report(path, str(exc)) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="只读检查 Harker/TAS 主量元素输入，并输出 JSON 报告。" + ) + parser.add_argument( + "input", type=Path, help=".csv、.txt 或 .xlsx 输入文件" + ) + parser.add_argument( + "--sheet", help="Excel 工作表名称，或从 0 开始的编号" + ) + parser.add_argument("--sample-column", help="明确指定样品编号列") + parser.add_argument("--group-column", help="明确指定可选分组列") + return parser.parse_args() + + +def main() -> int: + if hasattr(sys.stdout, "reconfigure"): + sys.stdout.reconfigure(encoding="utf-8") + if hasattr(sys.stderr, "reconfigure"): + sys.stderr.reconfigure(encoding="utf-8") + args = parse_args() + report = inspect_major_path( + args.input, + args.sheet, + args.sample_column, + args.group_column, + ) + json.dump(report, sys.stdout, ensure_ascii=False, indent=2) + sys.stdout.write("\n") + print(f"主量元素数据检查完成：{report['status']}", file=sys.stderr) + if report["status"] == "ready": + return 0 + if report["status"] == "error": + return 1 + return 2 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/skills/geoskills/scripts/plot_geochem_common.py b/skills/geoskills/scripts/plot_geochem_common.py new file mode 100644 index 0000000..660bd84 --- /dev/null +++ b/skills/geoskills/scripts/plot_geochem_common.py @@ -0,0 +1,239 @@ +#!/usr/bin/env python3 +"""Shared deterministic helpers for publication geochemical scatter figures.""" + +from __future__ import annotations + +import hashlib +from pathlib import Path +from typing import Any + +import numpy as np +import pandas as pd + +from plot_ree import FORMATS, GROUP_COLORS, MARKERS, PlottingError + + +SUBSCRIPT_TRANSLATION = str.maketrans("0123456789", "₀₁₂₃₄₅₆₇₈₉") + + +def analyte_label(analyte: str, unit: str) -> str: + """Return a compact final-size axis label with Unicode subscripts.""" + formula = analyte.translate(SUBSCRIPT_TRANSLATION) + return f"{formula} ({unit})" + + +def clean_linear_limits( + values: np.ndarray | pd.Series, + margin_fraction: float = 0.06, + nonnegative: bool = True, +) -> dict[str, float | str]: + """Return clean outward-rounded linear limits without clipping data.""" + if not 0.0 <= margin_fraction <= 0.25: + raise PlottingError("线性坐标边距必须在 0–0.25 之间。") + numeric = np.asarray(values, dtype=float).ravel() + valid = numeric[np.isfinite(numeric)] + if valid.size == 0: + raise PlottingError("没有可用于坐标范围计算的有限数值。") + data_min = float(valid.min()) + data_max = float(valid.max()) + if nonnegative and data_min < 0: + raise PlottingError("检测到负浓度，不能安全生成地球化学散点图。") + + raw_span = data_max - data_min + reference = max(abs(data_min), abs(data_max), 1.0) + span = raw_span if raw_span > 0 else reference * 0.1 + padded_min = data_min - span * margin_fraction + padded_max = data_max + span * margin_fraction + target_step = max((padded_max - padded_min) / 5.0, np.finfo(float).eps) + exponent = float(np.floor(np.log10(target_step))) + fraction = target_step / (10.0**exponent) + if fraction <= 1.0: + nice_fraction = 1.0 + elif fraction <= 2.0: + nice_fraction = 2.0 + elif fraction <= 2.5: + nice_fraction = 2.5 + elif fraction <= 5.0: + nice_fraction = 5.0 + else: + nice_fraction = 10.0 + step = float(nice_fraction * (10.0**exponent)) + + lower = float(np.floor(padded_min / step) * step) + upper = float(np.ceil(padded_max / step) * step) + if nonnegative and lower < 0: + lower = 0.0 + if np.isclose(lower, upper): + upper = lower + step + if lower > data_min or upper < data_max: + raise PlottingError("整洁坐标范围会裁切数据，已停止绘图。") + if np.isclose(lower, 0.0): + lower = 0.0 + if np.isclose(upper, 0.0): + upper = 0.0 + return { + "lower": lower, + "upper": upper, + "step": step, + "data_min": data_min, + "data_max": data_max, + "margin_fraction": margin_fraction, + "policy": "adaptive_linear_clean_outward_bounds", + } + + +def resolve_analyte_selection( + value: str, + available: list[str], + option_name: str, +) -> str: + """Resolve one case-insensitive canonical analyte selection.""" + canonical = {item.casefold(): item for item in available} + selected = canonical.get(value.strip().casefold()) + if selected is None: + raise PlottingError( + f"{option_name} 变量未在已验证输入中出现：{value}。" + ) + return selected + + +def resolve_analyte_list( + value: str | None, + available: list[str], + option_name: str, +) -> list[str]: + """Resolve a comma-separated case-insensitive analyte list.""" + if value is None: + return [] + requested = [item.strip() for item in value.split(",") if item.strip()] + if not requested: + raise PlottingError(f"{option_name} 不能为空。") + canonical = {item.casefold(): item for item in available} + selected: list[str] = [] + for item in requested: + resolved = canonical.get(item.casefold()) + if resolved is None: + raise PlottingError( + f"{option_name} 变量未在已验证输入中出现：{item}。" + ) + selected.append(resolved) + if len(set(selected)) != len(selected): + raise PlottingError(f"{option_name} 中不能重复同一变量。") + return selected + + +def filter_requested_groups( + frame: pd.DataFrame, + group_column: str | None, + selection: str | None, +) -> tuple[pd.DataFrame, list[str]]: + """Filter explicit groups without silently guessing category names.""" + if group_column is None: + if selection: + raise PlottingError("输入没有已验证分组列，不能使用 --groups。") + return frame.copy(), [] + working = frame.copy() + group_text = ( + working[group_column] + .astype("string") + .str.strip() + .fillna("Unspecified") + .replace("", "Unspecified") + ) + working[group_column] = group_text + available = [ + str(value) + for value in group_text.dropna().drop_duplicates().tolist() + if str(value) + ] + if selection is None: + selected = available + else: + requested = [ + item.strip() for item in selection.split(",") if item.strip() + ] + missing = [item for item in requested if item not in available] + if missing: + raise PlottingError( + "所选分组未在输入中出现：" + ", ".join(missing) + "。" + ) + selected = requested + filtered = working.loc[group_text.isin(selected)].copy() + if filtered.empty: + raise PlottingError("分组筛选后没有可绘制样品。") + return filtered, selected + + +def group_style_map(groups: list[str]) -> dict[str, dict[str, Any]]: + """Return deterministic colour-plus-marker styles for groups.""" + return { + group: { + "color": GROUP_COLORS[index % len(GROUP_COLORS)], + "marker": MARKERS[index % len(MARKERS)], + } + for index, group in enumerate(groups) + } + + +def validate_export_parameters( + width_mm: float, + height_mm: float, + dpi: int, +) -> None: + """Validate final-size and raster export parameters.""" + if not 50 <= width_mm <= 400 or not 50 <= height_mm <= 400: + raise PlottingError("图宽和图高必须在 50–400 mm 之间。") + if not 72 <= dpi <= 1200: + raise PlottingError("PNG/TIFF 分辨率必须在 72–1200 dpi 之间。") + + +def output_targets( + output_dir: Path, + stem: str, +) -> tuple[list[Path], Path, Path]: + """Resolve a complete figure bundle without accepting path-like stems.""" + if not stem or Path(stem).name != stem or Path(stem).suffix: + raise PlottingError("输出名称必须是不含路径和扩展名的文件名。") + figures = [ + output_dir / f"{stem}.{extension}" for extension in FORMATS + ] + source_data = output_dir / f"{stem}.source_data.csv" + report = output_dir / f"{stem}.report.json" + return figures, source_data, report + + +def ensure_outputs_available( + targets: list[Path], + overwrite: bool, +) -> None: + """Require an explicit overwrite decision for a complete bundle.""" + existing = [path for path in targets if path.exists()] + if existing and not overwrite: + raise PlottingError( + "输出文件已经存在；如需替换，请显式使用 --overwrite。" + ) + + +def save_figure_bundle( + figure: Any, + paths: list[Path], + dpi: int, +) -> None: + """Save vector and raster outputs from the same Matplotlib figure.""" + for path in paths: + options: dict[str, Any] = {"facecolor": "white"} + if path.suffix.lower() in {".png", ".tiff"}: + options["dpi"] = dpi + if path.suffix.lower() == ".tiff": + options["pil_kwargs"] = {"compression": "tiff_lzw"} + figure.savefig(path, **options) + + +def shareable_file_record(path: Path) -> dict[str, Any]: + """Describe an output without exposing its local directory.""" + return { + "format": path.suffix.lower().lstrip("."), + "filename": path.name, + "bytes": path.stat().st_size, + "sha256": hashlib.sha256(path.read_bytes()).hexdigest(), + } diff --git a/skills/geoskills/scripts/plot_harker.py b/skills/geoskills/scripts/plot_harker.py new file mode 100644 index 0000000..46989e7 --- /dev/null +++ b/skills/geoskills/scripts/plot_harker.py @@ -0,0 +1,770 @@ +#!/usr/bin/env python3 +"""Create publication-oriented Harker variation-diagram grids.""" + +from __future__ import annotations + +import argparse +import json +import math +import string +import sys +import textwrap +from pathlib import Path +from typing import Any + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from matplotlib.lines import Line2D +from matplotlib.ticker import MultipleLocator + +from inspect_data import InspectionError, issue +from inspect_major_data import ( + inspect_major_frame, + prepare_geochem_frame, + read_major_table, +) +from plot_geochem_common import ( + analyte_label, + clean_linear_limits, + ensure_outputs_available, + filter_requested_groups, + group_style_map, + output_targets, + resolve_analyte_list, + resolve_analyte_selection, + save_figure_bundle, + shareable_file_record, + validate_export_parameters, +) +from plot_ree import GROUP_COLORS, MARKERS, PlottingError + + +plt.rcParams["font.sans-serif"] = [ + "DejaVu Sans", + "Arial", + "Liberation Sans", +] + +DEFAULT_X = "SiO2" +DEFAULT_Y_CANDIDATES = [ + "TiO2", + "Al2O3", + "Fe2O3T", + "FeOT", + "MgO", + "CaO", + "Na2O", + "K2O", + "P2O5", +] +MAX_PANELS = 9 +AXIS_LABEL_SIZE = 6.2 +TICK_LABEL_SIZE = 5.6 +PANEL_LABEL_SIZE = 7.0 +HARKER_LITHOLOGY_ALERT_TERMS = ( + "carbonatite", + "kimberlite", + "lamproite", +) + + +def default_harker_y(available: list[str], x_analyte: str) -> list[str]: + """Choose a conventional compact major-oxide grid.""" + selected: list[str] = [] + total_iron_selected = False + for analyte in DEFAULT_Y_CANDIDATES: + if analyte not in available or analyte == x_analyte: + continue + if analyte in {"Fe2O3T", "FeOT"}: + if total_iron_selected: + continue + total_iron_selected = True + selected.append(analyte) + return selected[:MAX_PANELS] + + +def automatic_height(panel_count: int, columns: int) -> float: + """Return a journal-compatible height for the panel grid.""" + rows = math.ceil(panel_count / columns) + return {1: 86.0, 2: 120.0, 3: 155.0}.get(rows, 170.0) + + +def automatic_columns(panel_count: int) -> int: + """Choose a compact grid without avoidable empty panels.""" + if panel_count == 9: + return 3 + if panel_count >= 7: + return 4 + if panel_count >= 5: + return 3 + if panel_count == 4: + return 2 + return panel_count + + +def potentially_mixed_harker_groups(groups: list[str]) -> list[str]: + """Flag specialist lithologies mixed with other Harker groups.""" + flagged = [ + group + for group in groups + if any( + term in group.casefold() + for term in HARKER_LITHOLOGY_ALERT_TERMS + ) + ] + if flagged and len(flagged) < len(groups): + return flagged + return [] + + +def add_shared_group_legend( + figure: Any, + handles: list[Any], + labels: list[str], + has_title: bool, +) -> dict[str, Any]: + """Fit a shared legend inside the figure canvas without clipping.""" + display_labels = [ + textwrap.fill( + label, + width=38, + break_long_words=False, + break_on_hyphens=False, + ) + for label in labels + ] + anchor_y = 0.945 if has_title else 0.985 + legend = None + legend_box = None + selected_columns = 1 + average_label_length = sum(len(label) for label in labels) / len(labels) + max_columns = min(4, len(labels)) + if len(labels) >= 4 and average_label_length > 24: + max_columns = min(2, len(labels)) + elif len(labels) == 3 and average_label_length > 30: + max_columns = 2 + for candidate_columns in range(max_columns, 0, -1): + legend = figure.legend( + handles, + display_labels, + loc="upper center", + bbox_to_anchor=(0.5, anchor_y), + ncol=candidate_columns, + frameon=False, + fontsize=5.8, + handletextpad=0.35, + columnspacing=0.9, + ) + figure.canvas.draw() + renderer = figure.canvas.get_renderer() + legend_box = legend.get_window_extent(renderer).transformed( + figure.transFigure.inverted() + ) + if legend_box.x0 >= 0.02 and legend_box.x1 <= 0.98: + selected_columns = candidate_columns + break + legend.remove() + legend = None + + if legend is None or legend_box is None: + raise PlottingError("共享图例无法安全放入图幅。") + axes_top = max(0.58, min(0.94, float(legend_box.y0) - 0.025)) + return { + "position": "shared_figure_top", + "columns": selected_columns, + "rows": math.ceil(len(labels) / selected_columns), + "axes_top": axes_top, + "within_figure": bool( + legend_box.x0 >= 0.0 + and legend_box.x1 <= 1.0 + and legend_box.y0 >= 0.0 + and legend_box.y1 <= 1.0 + ), + } + + +def build_harker_figure( + frame: pd.DataFrame, + sample_column: str, + group_column: str | None, + groups: list[str], + x_analyte: str, + y_analytes: list[str], + units: dict[str, str], + title: str | None, + width_mm: float, + height_mm: float, + columns: int, + axes_frame: str, + margin_fraction: float, +) -> tuple[Any, dict[str, Any]]: + """Build one Harker grid from validated canonical data.""" + x_limits = clean_linear_limits( + frame[x_analyte], margin_fraction=margin_fraction + ) + rows = math.ceil(len(y_analytes) / columns) + figure, axes = plt.subplots( + rows, + columns, + figsize=(width_mm / 25.4, height_mm / 25.4), + sharex=True, + squeeze=False, + ) + styles = group_style_map(groups) + panel_reports: list[dict[str, Any]] = [] + missing_pairs: dict[str, int] = {} + + for index, analyte in enumerate(y_analytes): + ax = axes.flat[index] + valid_pair_count = int( + ( + pd.to_numeric(frame[x_analyte], errors="coerce").notna() + & pd.to_numeric(frame[analyte], errors="coerce").notna() + ).sum() + ) + if valid_pair_count < 2: + raise PlottingError( + f"{x_analyte}–{analyte} 只有 {valid_pair_count} 个完整数据点；" + "至少需要 2 个。" + ) + y_limits = clean_linear_limits( + frame[analyte], margin_fraction=margin_fraction + ) + missing_pairs[analyte] = int(len(frame) - valid_pair_count) + + if group_column is None: + x_values = pd.to_numeric(frame[x_analyte], errors="coerce") + y_values = pd.to_numeric(frame[analyte], errors="coerce") + finite = x_values.notna() & y_values.notna() + ax.scatter( + x_values[finite], + y_values[finite], + s=22, + c="#3569A8", + marker="o", + edgecolors="white", + linewidths=0.45, + alpha=0.92, + zorder=3, + ) + else: + for group in groups: + subset = frame.loc[ + frame[group_column].astype(str) == group + ] + x_values = pd.to_numeric( + subset[x_analyte], errors="coerce" + ) + y_values = pd.to_numeric( + subset[analyte], errors="coerce" + ) + finite = x_values.notna() & y_values.notna() + style = styles[group] + ax.scatter( + x_values[finite], + y_values[finite], + s=22, + c=style["color"], + marker=style["marker"], + edgecolors="white", + linewidths=0.45, + alpha=0.92, + zorder=3, + ) + + ax.set_xlim( + float(x_limits["lower"]), float(x_limits["upper"]) + ) + ax.set_ylim( + float(y_limits["lower"]), float(y_limits["upper"]) + ) + ax.xaxis.set_major_locator( + MultipleLocator(float(x_limits["step"])) + ) + ax.yaxis.set_major_locator( + MultipleLocator(float(y_limits["step"])) + ) + ax.set_ylabel( + analyte_label(analyte, units[analyte]), + fontsize=AXIS_LABEL_SIZE, + labelpad=2.2, + ) + if index // columns != rows - 1: + ax.tick_params(axis="x", labelbottom=False) + ax.tick_params( + axis="both", + which="major", + direction="out", + length=2.8, + labelsize=TICK_LABEL_SIZE, + pad=1.8, + ) + ax.grid(False) + for spine in ax.spines.values(): + spine.set_linewidth(0.7) + spine.set_color("#1A1A1A") + show_full = axes_frame == "full" + ax.spines["top"].set_visible(show_full) + ax.spines["right"].set_visible(show_full) + ax.text( + 0.025, + 0.975, + string.ascii_lowercase[index], + transform=ax.transAxes, + ha="left", + va="top", + fontsize=PANEL_LABEL_SIZE, + fontweight="bold", + ) + panel_reports.append( + { + "panel": string.ascii_lowercase[index], + "x": x_analyte, + "y": analyte, + "complete_pairs": valid_pair_count, + "x_limits": x_limits, + "y_limits": y_limits, + } + ) + + for index in range(len(y_analytes), rows * columns): + axes.flat[index].set_axis_off() + + legend_info = { + "position": "none", + "columns": 0, + "rows": 0, + "axes_top": 0.90 if title else 0.94, + "within_figure": True, + } + if group_column is not None: + handles = [ + Line2D( + [], + [], + linestyle="none", + marker=styles[group]["marker"], + markerfacecolor=styles[group]["color"], + markeredgecolor="white", + markeredgewidth=0.45, + markersize=5.2, + ) + for group in groups + ] + labels = [ + f"{group} (n={int((frame[group_column].astype(str) == group).sum())})" + for group in groups + ] + legend_info = add_shared_group_legend( + figure, + handles, + labels, + has_title=bool(title), + ) + + figure.subplots_adjust( + left=0.075, + right=0.985, + bottom=0.11, + top=legend_info["axes_top"], + hspace=0.34, + wspace=0.46 if columns == 4 else 0.34, + ) + figure.supxlabel( + analyte_label(x_analyte, units[x_analyte]), + x=0.53, + y=0.025, + fontsize=AXIS_LABEL_SIZE, + ) + if title: + figure.suptitle(title, y=0.995, fontsize=8) + return figure, { + "panel_count": len(y_analytes), + "rows": rows, + "columns": columns, + "sample_count": int(len(frame)), + "group_count": len(groups) if group_column is not None else 0, + "legend_position": legend_info["position"], + "legend_columns": legend_info["columns"], + "legend_rows": legend_info["rows"], + "legend_within_figure": legend_info["within_figure"], + "shared_x_label": True, + "axes_frame": axes_frame, + "x_limits": x_limits, + "panels": panel_reports, + "missing_pairs": missing_pairs, + "palette_repeated": len(groups) > len(GROUP_COLORS), + "marker_repeated": len(groups) > len(MARKERS), + "colour_is_not_the_only_identifier": True, + "group_encoding": ( + "colour plus marker shape" + if group_column is not None + else "not applicable" + ), + } + + +def harker_error(path: Path, message: str) -> dict[str, Any]: + return { + "status": "error", + "operation": "harker_plot", + "source": {"format": path.suffix.lower()}, + "issues": [issue("E650", "error", message)], + } + + +def plot_harker_path( + input_path: Path, + output_dir: Path, + stem: str | None = None, + requested_sheet: str | None = None, + requested_sample_column: str | None = None, + requested_group_column: str | None = None, + requested_x: str = DEFAULT_X, + requested_y: str | None = None, + requested_groups: str | None = None, + title: str | None = None, + width_mm: float = 183.0, + height_mm: float | None = None, + columns: int | None = None, + dpi: int = 600, + axes_frame: str = "full", + margin_fraction: float = 0.06, + overwrite: bool = False, +) -> dict[str, Any]: + """Validate input and export a Harker figure bundle.""" + figure = None + try: + if columns is not None and not 1 <= columns <= 4: + raise PlottingError("--columns 必须在 1–4 之间。") + frame, source = read_major_table(input_path, requested_sheet) + if frame is None: + return { + "status": "needs_sheet", + "operation": "harker_plot", + "source": source, + "issues": [ + issue( + "E111", + "review", + "Excel 文件包含多个工作表，请使用 --sheet 指定一个工作表。", + sheet_names=source["sheet_names"], + ) + ], + } + inspection = inspect_major_frame( + frame, + source, + requested_sample_column, + requested_group_column, + ) + if inspection["status"] != "ready": + return { + "status": "blocked", + "operation": "harker_plot", + "input_inspection": inspection, + "issues": [ + issue( + "E651", + "review", + "输入数据未通过安全检查，因此没有生成图像。", + ) + ], + } + + available = [ + item["analyte"] + for item in inspection["analytes"]["recognized"] + ] + units = { + item["analyte"]: item["unit"] + for item in inspection["analytes"]["recognized"] + } + x_analyte = resolve_analyte_selection( + requested_x, available, "--x" + ) + y_analytes = resolve_analyte_list( + requested_y, available, "--y" + ) + if not y_analytes: + y_analytes = default_harker_y(available, x_analyte) + if not y_analytes: + raise PlottingError( + "没有可自动选择的 Harker 纵轴变量，请使用 --y 明确指定。" + ) + if x_analyte in y_analytes: + raise PlottingError("Harker 图的 X 和 Y 变量不能相同。") + if len(y_analytes) > MAX_PANELS: + raise PlottingError( + f"单张 Harker 图最多支持 {MAX_PANELS} 个纵轴面板。" + ) + resolved_columns = ( + columns + if columns is not None + else automatic_columns(len(y_analytes)) + ) + + canonical = prepare_geochem_frame(frame, inspection) + sample_column = inspection["sample_id_candidates"][0] + group_column = ( + inspection["group_candidates"][0] + if len(inspection["group_candidates"]) == 1 + else None + ) + canonical, groups = filter_requested_groups( + canonical, group_column, requested_groups + ) + resolved_height = ( + height_mm + if height_mm is not None + else automatic_height(len(y_analytes), resolved_columns) + ) + validate_export_parameters(width_mm, resolved_height, dpi) + resolved_stem = stem or f"harker_{x_analyte.lower()}" + figure_paths, source_path, report_path = output_targets( + output_dir, resolved_stem + ) + ensure_outputs_available( + [*figure_paths, source_path, report_path], overwrite + ) + + figure, plot_info = build_harker_figure( + canonical, + sample_column, + group_column, + groups, + x_analyte, + y_analytes, + units, + title, + width_mm, + resolved_height, + resolved_columns, + axes_frame, + margin_fraction, + ) + output_dir.mkdir(parents=True, exist_ok=True) + save_figure_bundle(figure, figure_paths, dpi) + source_columns = [sample_column] + if group_column is not None: + source_columns.append(group_column) + source_columns.extend([x_analyte, *y_analytes]) + canonical.loc[:, source_columns].to_csv( + source_path, + index=False, + encoding="utf-8", + float_format="%.10g", + ) + plt.close(figure) + figure = None + + run_issues = list(inspection["issues"]) + incomplete = { + analyte: count + for analyte, count in plot_info["missing_pairs"].items() + if count + } + if incomplete: + run_issues.append( + issue( + "W651", + "warning", + "部分 Harker 面板存在缺失配对，已跳过相应散点。", + missing_pairs=incomplete, + ) + ) + if plot_info["palette_repeated"] or plot_info["marker_repeated"]: + run_issues.append( + issue( + "W652", + "warning", + "分组数量超过基础颜色或符号数量；建议筛选分组或分图。", + group_count=plot_info["group_count"], + ) + ) + mixed_groups = potentially_mixed_harker_groups(groups) + if mixed_groups: + run_issues.append( + issue( + "W653", + "warning", + "分组名称提示图中混合了差异很大的岩石类型；" + "解释成岩浆演化趋势前应筛选成分与成因上可比较的样品组。", + flagged_groups=mixed_groups, + ) + ) + + report = { + "status": "ready", + "operation": "harker_plot", + "source": source, + "figure_contract": { + "core_conclusion": ( + "Show how selected whole-rock analytes covary with one " + "declared differentiation index across samples and groups." + ), + "archetype": "quantitative grid", + "backend": "Python/matplotlib", + "role": "comparative and discovery evidence", + "evidence": "paired whole-rock concentrations", + "target_output": "double-column publication figure", + "review_risks": [ + "unit ambiguity", + "missing paired values", + "mixed lithologies or unrelated magma series", + "correlation interpreted as process", + "overplotting and group confusion", + ], + }, + "configuration": { + "sample_column": sample_column, + "group_column": group_column, + "groups": groups, + "x": x_analyte, + "y": y_analytes, + "units": {key: units[key] for key in [x_analyte, *y_analytes]}, + "trend_lines": "none", + "axes_frame": axes_frame, + "linear_limit_policy": ( + "adaptive_linear_clean_outward_bounds" + ), + "margin_fraction": margin_fraction, + "columns": resolved_columns, + "columns_mode": ( + "explicit" if columns is not None else "automatic" + ), + "width_mm": width_mm, + "height_mm": resolved_height, + "png_dpi": dpi, + "tiff_dpi": dpi, + }, + "plot": plot_info, + "outputs": [ + shareable_file_record(path) for path in figure_paths + ], + "source_data": shareable_file_record(source_path), + "submission_qa": { + "final_size_mm": [width_mm, resolved_height], + "svg_text_editable": True, + "pdf_font_type": 42, + "raster_dpi": dpi, + "tiff_compression": "LZW", + "white_background": True, + "shared_legend": group_column is not None, + "colourblind_support": "group colour plus marker shape", + "source_data_exported": True, + }, + "report_file": report_path.name, + "issues": run_issues, + "interpretation_guidance": [ + "Describe direction, curvature, scatter, clusters, and outliers before proposing a process.", + "Check whether compared samples belong to a coherent suite before interpreting a differentiation trend.", + "Consider alteration, analytical uncertainty, mixing, assimilation, and source variation as alternatives.", + "Do not treat an apparent linear relation as proof of fractional crystallization.", + ], + "scientific_caveat": ( + "Harker covariance alone does not identify a unique " + "petrogenetic process or establish genetic relationships." + ), + } + report_path.write_text( + json.dumps(report, ensure_ascii=False, indent=2) + "\n", + encoding="utf-8", + ) + return report + except ( + InspectionError, + PlottingError, + OSError, + ValueError, + ) as exc: + return harker_error(input_path, str(exc)) + finally: + if figure is not None: + plt.close(figure) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="生成可自定义 X/Y 变量的投稿级 Harker 图解。" + ) + parser.add_argument("input", type=Path, help="CSV、TXT 或 Excel 输入表格") + parser.add_argument("--output-dir", type=Path, required=True) + parser.add_argument("--stem", help="不含扩展名的输出文件名") + parser.add_argument("--sheet", help="Excel 工作表名称或编号") + parser.add_argument("--sample-column", help="明确指定样品编号列") + parser.add_argument("--group-column", help="明确指定可选分组列") + parser.add_argument("--x", default=DEFAULT_X, help="横轴变量，默认 SiO2") + parser.add_argument( + "--y", help="逗号分隔的纵轴变量；不指定时自动选择常用主量氧化物" + ) + parser.add_argument("--groups", help="逗号分隔的待绘制分组") + parser.add_argument("--title", help="可选图题") + parser.add_argument("--width-mm", type=float, default=183.0) + parser.add_argument( + "--height-mm", + type=float, + help="图高；不指定时按面板数自动选择投稿尺寸", + ) + parser.add_argument( + "--columns", + type=int, + help="每行面板数（1–4）；不指定时自动选择紧凑布局", + ) + parser.add_argument("--dpi", type=int, default=600) + parser.add_argument( + "--axes-frame", + choices=("open", "full"), + default="full", + ) + parser.add_argument( + "--margin", + type=float, + default=0.06, + help="线性坐标数据边距比例，默认 0.06", + ) + parser.add_argument("--overwrite", action="store_true") + return parser.parse_args() + + +def main() -> int: + if hasattr(sys.stdout, "reconfigure"): + sys.stdout.reconfigure(encoding="utf-8") + if hasattr(sys.stderr, "reconfigure"): + sys.stderr.reconfigure(encoding="utf-8") + args = parse_args() + report = plot_harker_path( + args.input, + args.output_dir, + stem=args.stem, + requested_sheet=args.sheet, + requested_sample_column=args.sample_column, + requested_group_column=args.group_column, + requested_x=args.x, + requested_y=args.y, + requested_groups=args.groups, + title=args.title, + width_mm=args.width_mm, + height_mm=args.height_mm, + columns=args.columns, + dpi=args.dpi, + axes_frame=args.axes_frame, + margin_fraction=args.margin, + overwrite=args.overwrite, + ) + json.dump(report, sys.stdout, ensure_ascii=False, indent=2) + sys.stdout.write("\n") + print(f"Harker 绘图完成：{report['status']}", file=sys.stderr) + if report["status"] == "ready": + return 0 + if report["status"] == "error": + return 1 + return 2 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/skills/geoskills/scripts/plot_tas.py b/skills/geoskills/scripts/plot_tas.py new file mode 100644 index 0000000..55e5e7c --- /dev/null +++ b/skills/geoskills/scripts/plot_tas.py @@ -0,0 +1,863 @@ +#!/usr/bin/env python3 +"""Create and classify publication-oriented volcanic TAS diagrams.""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import sys +from collections import Counter +from pathlib import Path +from typing import Any + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from matplotlib.lines import Line2D +from matplotlib.patches import Polygon +from matplotlib.path import Path as MatplotlibPath + +from inspect_data import InspectionError, issue +from inspect_major_data import ( + inspect_major_frame, + prepare_geochem_frame, + read_major_table, +) +from plot_geochem_common import ( + ensure_outputs_available, + filter_requested_groups, + group_style_map, + output_targets, + save_figure_bundle, + shareable_file_record, + validate_export_parameters, +) +from plot_ree import GROUP_COLORS, MARKERS, PlottingError, configure_boxed_legend + + +plt.rcParams["font.sans-serif"] = [ + "DejaVu Sans", + "Arial", + "Liberation Sans", +] + +SKILL_DIR = Path(__file__).resolve().parents[1] +DEFAULT_MODEL_PATH = ( + SKILL_DIR / "assets" / "classification" / "tas-lemaitre-2002.json" +) +COMPOSITION_BASES = ("unknown", "anhydrous-normalized", "as-reported") +LEGEND_LAYOUTS = ("inside-auto", "outside") +POTENTIALLY_INCOMPATIBLE_GROUP_TERMS = ( + "syenite", + "granite", + "gabbro", + "diorite", + "carbonatite", + "plutonic", +) + + +def load_tas_model(path: Path) -> dict[str, Any]: + """Load and minimally validate the versioned TAS polygon model.""" + try: + model = json.loads(path.read_text(encoding="utf-8")) + except (OSError, json.JSONDecodeError) as exc: + raise PlottingError(f"无法读取 TAS 边界资产：{exc}") from exc + if model.get("diagram_type") != "total_alkali_silica": + raise PlottingError("TAS 边界资产类型不正确。") + fields = model.get("fields") + if not isinstance(fields, list) or len(fields) != 15: + raise PlottingError("TAS 边界资产必须包含 15 个火山岩字段。") + ids = [field.get("id") for field in fields] + if len(set(ids)) != len(ids): + raise PlottingError("TAS 边界资产包含重复字段 ID。") + for field in fields: + polygon = field.get("polygon") + if not isinstance(polygon, list) or len(polygon) < 3: + raise PlottingError( + f"TAS 字段 {field.get('id')} 缺少有效多边形。" + ) + return model + + +def tas_model_summary(path: Path, model: dict[str, Any]) -> dict[str, Any]: + """Return shareable provenance for the TAS boundary asset.""" + return { + "id": model["id"], + "display_name": model["display_name"], + "classification_doi": model["source"]["classification_doi"], + "boundary_doi": model["source"]["boundary_doi"], + "asset": path.name, + "sha256": hashlib.sha256(path.read_bytes()).hexdigest(), + "crosscheck": model.get("crosscheck"), + "boundary_policy": model["boundary_policy"], + } + + +def point_segment_distance( + point: np.ndarray, + start: np.ndarray, + end: np.ndarray, +) -> float: + """Return Euclidean distance from a point to one line segment.""" + vector = end - start + length_squared = float(np.dot(vector, vector)) + if length_squared == 0: + return float(np.linalg.norm(point - start)) + fraction = float(np.dot(point - start, vector) / length_squared) + fraction = min(1.0, max(0.0, fraction)) + projection = start + fraction * vector + return float(np.linalg.norm(point - projection)) + + +def point_on_boundary( + x: float, + y: float, + polygon: list[list[float]], + tolerance: float, +) -> bool: + """Return whether a point lies on any polygon segment.""" + point = np.asarray([x, y], dtype=float) + vertices = np.asarray(polygon, dtype=float) + for index in range(len(vertices)): + start = vertices[index] + end = vertices[(index + 1) % len(vertices)] + if point_segment_distance(point, start, end) <= tolerance: + return True + return False + + +def classify_tas_point( + x: float, + y: float, + model: dict[str, Any], +) -> dict[str, str | None]: + """Classify one point, preserving boundary and outside states.""" + x_limits = model["axes"]["x"]["limits"] + y_limits = model["axes"]["y"]["limits"] + if not np.isfinite(x) or not np.isfinite(y): + return { + "status": "missing", + "field": None, + "name": None, + } + if ( + x < x_limits[0] + or x > x_limits[1] + or y < y_limits[0] + or y > y_limits[1] + ): + return { + "status": "outside_model", + "field": None, + "name": None, + } + tolerance = float( + model["boundary_policy"]["tolerance_wt_percent"] + ) + boundary_fields = [ + field["id"] + for field in model["fields"] + if point_on_boundary(x, y, field["polygon"], tolerance) + ] + if boundary_fields: + return { + "status": "boundary_review", + "field": "/".join(boundary_fields), + "name": None, + } + + matches: list[dict[str, Any]] = [] + for field in model["fields"]: + vertices = np.asarray(field["polygon"], dtype=float) + closed = np.vstack([vertices, vertices[0]]) + if MatplotlibPath(closed).contains_point((x, y)): + matches.append(field) + if len(matches) == 1: + return { + "status": "classified", + "field": matches[0]["id"], + "name": matches[0]["name"], + } + if len(matches) > 1: + return { + "status": "ambiguous_review", + "field": "/".join(field["id"] for field in matches), + "name": None, + } + return { + "status": "outside_fields", + "field": None, + "name": None, + } + + +def classify_tas_frame( + frame: pd.DataFrame, + sample_column: str, + group_column: str | None, + model: dict[str, Any], +) -> pd.DataFrame: + """Return TAS coordinates and deterministic classification states.""" + output = frame.loc[:, [sample_column]].copy() + if group_column is not None: + output[group_column] = frame[group_column] + output["SiO2"] = pd.to_numeric(frame["SiO2"], errors="coerce") + output["Na2O"] = pd.to_numeric(frame["Na2O"], errors="coerce") + output["K2O"] = pd.to_numeric(frame["K2O"], errors="coerce") + output["TotalAlkali"] = output["Na2O"] + output["K2O"] + classifications = [ + classify_tas_point(float(x), float(y), model) + for x, y in zip(output["SiO2"], output["TotalAlkali"]) + ] + output["TAS_status"] = [item["status"] for item in classifications] + output["TAS_field"] = [item["field"] for item in classifications] + output["TAS_name"] = [item["name"] for item in classifications] + return output + + +def legend_overlaps_content( + ax: Any, + legend: Any, + points: np.ndarray, + field_texts: list[Any], +) -> bool: + """Return whether a legend overlaps sample points or field labels.""" + ax.figure.canvas.draw() + renderer = ax.figure.canvas.get_renderer() + box = legend.get_window_extent(renderer).expanded(1.04, 1.08) + if points.size: + display = ax.transData.transform(points) + inside = ( + (display[:, 0] >= box.x0) + & (display[:, 0] <= box.x1) + & (display[:, 1] >= box.y0) + & (display[:, 1] <= box.y1) + ) + if inside.any(): + return True + return any( + box.overlaps(text.get_window_extent(renderer).expanded(1.02, 1.06)) + for text in field_texts + ) + + +def add_tas_legend( + figure: Any, + ax: Any, + handles: list[Any], + labels: list[str], + points: np.ndarray, + field_texts: list[Any], + requested_layout: str, +) -> tuple[Any, str, bool]: + """Place a collision-checked TAS legend or safely move it outside.""" + if requested_layout == "inside-auto": + candidates = [ + "lower right", + "upper right", + "lower left", + "upper left", + ] + for location in candidates: + legend = ax.legend( + handles, + labels, + loc=location, + fontsize=5.7, + handletextpad=0.35, + labelspacing=0.35, + borderpad=0.45, + frameon=True, + fancybox=False, + ) + configure_boxed_legend(legend) + if not legend_overlaps_content( + ax, legend, points, field_texts + ): + return legend, f"inside_{location.replace(' ', '_')}", False + legend.remove() + figure.subplots_adjust(right=0.75) + legend = ax.legend( + handles, + labels, + loc="center left", + bbox_to_anchor=(1.02, 0.5), + fontsize=5.8, + handletextpad=0.4, + labelspacing=0.38, + frameon=False, + ) + return legend, "outside_right_fallback", requested_layout == "inside-auto" + + +def build_tas_figure( + classified: pd.DataFrame, + sample_column: str, + group_column: str | None, + groups: list[str], + model: dict[str, Any], + composition_basis: str, + title: str | None, + width_mm: float, + height_mm: float, + legend_layout: str, +) -> tuple[Any, dict[str, Any]]: + """Build one fixed-geometry volcanic TAS figure.""" + figure, ax = plt.subplots( + figsize=(width_mm / 25.4, height_mm / 25.4) + ) + figure.subplots_adjust( + left=0.105, + right=0.965, + bottom=0.15, + top=0.93 if title else 0.965, + ) + field_texts: list[Any] = [] + fills = ("#FAFAFA", "#F2F2F2") + for index, field in enumerate(model["fields"]): + patch = Polygon( + field["polygon"], + closed=True, + facecolor=fills[index % len(fills)], + edgecolor="#686868", + linewidth=0.55, + zorder=1, + ) + ax.add_patch(patch) + field_texts.append( + ax.text( + field["label_position"][0], + field["label_position"][1], + field["plot_label"], + ha="center", + va="center", + fontsize=5.0, + color="#4D4D4D", + linespacing=0.9, + zorder=2, + ) + ) + + styles = group_style_map(groups) + if group_column is None: + finite = ( + classified["SiO2"].notna() + & classified["TotalAlkali"].notna() + ) + ax.scatter( + classified.loc[finite, "SiO2"], + classified.loc[finite, "TotalAlkali"], + s=28, + c="#3569A8", + marker="o", + edgecolors="white", + linewidths=0.55, + alpha=0.94, + zorder=4, + ) + legend_position = "none" + legend_fallback = False + else: + for group in groups: + subset = classified.loc[ + classified[group_column].astype(str) == group + ] + finite = ( + subset["SiO2"].notna() + & subset["TotalAlkali"].notna() + ) + style = styles[group] + ax.scatter( + subset.loc[finite, "SiO2"], + subset.loc[finite, "TotalAlkali"], + s=28, + c=style["color"], + marker=style["marker"], + edgecolors="white", + linewidths=0.55, + alpha=0.94, + zorder=4, + ) + handles = [ + Line2D( + [], + [], + linestyle="none", + marker=styles[group]["marker"], + markerfacecolor=styles[group]["color"], + markeredgecolor="white", + markeredgewidth=0.5, + markersize=5.5, + ) + for group in groups + ] + labels = [ + f"{group} (n={int((classified[group_column].astype(str) == group).sum())})" + for group in groups + ] + finite_points = classified.loc[ + classified["SiO2"].notna() + & classified["TotalAlkali"].notna(), + ["SiO2", "TotalAlkali"], + ].to_numpy(dtype=float) + _, legend_position, legend_fallback = add_tas_legend( + figure, + ax, + handles, + labels, + finite_points, + field_texts, + legend_layout, + ) + + ax.set_xlim(model["axes"]["x"]["limits"]) + ax.set_ylim(model["axes"]["y"]["limits"]) + ax.set_xlabel("SiO₂ (wt%)") + ax.set_ylabel("Na₂O + K₂O (wt%)") + ax.set_xticks(np.arange(35, 91, 5)) + ax.set_yticks(np.arange(0, 21, 2)) + ax.tick_params( + axis="both", which="major", direction="out", length=3 + ) + ax.grid(False) + for spine in ax.spines.values(): + spine.set_visible(True) + spine.set_color("#1A1A1A") + spine.set_linewidth(0.75) + if title: + ax.set_title(title, fontsize=8, pad=6) + basis_label = ( + "declared anhydrous 100% basis" + if composition_basis == "anhydrous-normalized" + else "as-reported values; classification is provisional" + ) + figure.text( + 0.105, + 0.045, + "IUGS volcanic TAS (Le Maitre et al., 2002); " + + basis_label + + "; boundary points require review.", + fontsize=5.6, + color="#4D4D4D", + ha="left", + ) + counts = Counter(classified["TAS_status"].tolist()) + return figure, { + "sample_count": int(len(classified)), + "group_count": len(groups) if group_column is not None else 0, + "classification_status_counts": dict(counts), + "field_counts": { + str(key): int(value) + for key, value in classified["TAS_field"] + .dropna() + .value_counts() + .items() + }, + "legend_layout_requested": legend_layout, + "legend_position": legend_position, + "legend_fallback": legend_fallback, + "axes_frame": "full", + "x_limits": model["axes"]["x"]["limits"], + "y_limits": model["axes"]["y"]["limits"], + "palette_repeated": len(groups) > len(GROUP_COLORS), + "marker_repeated": len(groups) > len(MARKERS), + "colour_is_not_the_only_identifier": True, + "group_encoding": ( + "colour plus marker shape" + if group_column is not None + else "not applicable" + ), + } + + +def tas_error(path: Path, message: str) -> dict[str, Any]: + return { + "status": "error", + "operation": "tas_plot", + "source": {"format": path.suffix.lower()}, + "issues": [issue("E700", "error", message)], + } + + +def plot_tas_path( + input_path: Path, + output_dir: Path, + stem: str | None = None, + requested_sheet: str | None = None, + requested_sample_column: str | None = None, + requested_group_column: str | None = None, + requested_groups: str | None = None, + confirm_volcanic: bool = False, + composition_basis: str = "unknown", + title: str | None = None, + width_mm: float = 150.0, + height_mm: float = 120.0, + dpi: int = 600, + legend_layout: str = "inside-auto", + overwrite: bool = False, + model_path: Path = DEFAULT_MODEL_PATH, +) -> dict[str, Any]: + """Validate, classify, and export a volcanic TAS figure bundle.""" + figure = None + try: + if composition_basis not in COMPOSITION_BASES: + raise PlottingError( + "未知 composition basis；请使用 anhydrous-normalized 或 as-reported。" + ) + if legend_layout not in LEGEND_LAYOUTS: + raise PlottingError("未知 TAS 图例布局。") + validate_export_parameters(width_mm, height_mm, dpi) + model = load_tas_model(model_path) + frame, source = read_major_table(input_path, requested_sheet) + if frame is None: + return { + "status": "needs_sheet", + "operation": "tas_plot", + "source": source, + "model": tas_model_summary(model_path, model), + "issues": [ + issue( + "E111", + "review", + "Excel 文件包含多个工作表，请使用 --sheet 指定一个工作表。", + sheet_names=source["sheet_names"], + ) + ], + } + inspection = inspect_major_frame( + frame, + source, + requested_sample_column, + requested_group_column, + ) + gate_issues: list[dict[str, Any]] = [] + if not confirm_volcanic: + gate_issues.append( + issue( + "E711", + "review", + "TAS 仅用于本版本声明的火山岩分类；请确认样品为火山岩。", + ) + ) + if composition_basis == "unknown": + gate_issues.append( + issue( + "E712", + "review", + "请声明数据是 anhydrous-normalized 还是 as-reported。", + ) + ) + if inspection["status"] != "ready" or gate_issues: + return { + "status": "blocked", + "operation": "tas_plot", + "model": tas_model_summary(model_path, model), + "input_inspection": inspection, + "issues": [ + *gate_issues, + *( + [ + issue( + "E713", + "review", + "输入数据未通过安全检查，因此没有生成 TAS 图。", + ) + ] + if inspection["status"] != "ready" + else [] + ), + ], + } + + available = { + item["analyte"] + for item in inspection["analytes"]["recognized"] + } + required = {"SiO2", "Na2O", "K2O"} + missing = sorted(required - available) + if missing: + raise PlottingError( + "TAS 缺少必要主量氧化物：" + ", ".join(missing) + "。" + ) + units = { + item["analyte"]: item["unit"] + for item in inspection["analytes"]["recognized"] + } + if any(units[analyte] != "wt%" for analyte in required): + raise PlottingError("TAS 的 SiO2、Na2O 和 K2O 必须全部为 wt%。") + + canonical = prepare_geochem_frame(frame, inspection) + sample_column = inspection["sample_id_candidates"][0] + group_column = ( + inspection["group_candidates"][0] + if len(inspection["group_candidates"]) == 1 + else None + ) + canonical, groups = filter_requested_groups( + canonical, group_column, requested_groups + ) + classified = classify_tas_frame( + canonical, + sample_column, + group_column, + model, + ) + resolved_stem = stem or "tas_volcanic" + figure_paths, source_path, report_path = output_targets( + output_dir, resolved_stem + ) + ensure_outputs_available( + [*figure_paths, source_path, report_path], overwrite + ) + figure, plot_info = build_tas_figure( + classified, + sample_column, + group_column, + groups, + model, + composition_basis, + title, + width_mm, + height_mm, + legend_layout, + ) + output_dir.mkdir(parents=True, exist_ok=True) + save_figure_bundle(figure, figure_paths, dpi) + classified.to_csv( + source_path, + index=False, + encoding="utf-8", + float_format="%.10g", + ) + plt.close(figure) + figure = None + + run_issues = list(inspection["issues"]) + if composition_basis == "as-reported": + run_issues.append( + issue( + "W711", + "warning", + "使用 as-reported 数值绘制；分类为临时结果，投稿前应复核无挥发分 100% 归一化基础。", + ) + ) + for status, code, message in [ + ( + "boundary_review", + "W712", + "部分样品落在 TAS 字段边界上，未强制指定岩类。", + ), + ( + "outside_model", + "W713", + "部分样品超出 TAS 模型固定坐标范围。", + ), + ( + "outside_fields", + "W714", + "部分样品位于坐标范围内但不属于已定义字段。", + ), + ]: + count = plot_info["classification_status_counts"].get(status, 0) + if count: + run_issues.append( + issue(code, "warning", message, sample_count=count) + ) + if group_column is not None: + suspicious = [ + group + for group in groups + if any( + term in group.casefold() + for term in POTENTIALLY_INCOMPATIBLE_GROUP_TERMS + ) + ] + if suspicious: + run_issues.append( + issue( + "W715", + "warning", + "部分分组名称提示其可能不属于火山岩 TAS 适用范围，请人工复核。", + groups=suspicious, + ) + ) + if plot_info["legend_fallback"]: + run_issues.append( + issue( + "W716", + "warning", + "图内图例会遮挡样品点或分类标签，已改为右侧布局。", + ) + ) + if plot_info["palette_repeated"] or plot_info["marker_repeated"]: + run_issues.append( + issue( + "W717", + "warning", + "分组数量超过基础颜色或符号数量；建议筛选分组。", + group_count=plot_info["group_count"], + ) + ) + + report = { + "status": "ready", + "operation": "tas_plot", + "source": source, + "model": tas_model_summary(model_path, model), + "figure_contract": { + "core_conclusion": ( + "Place confirmed volcanic whole-rock analyses in declared " + "IUGS total-alkali–silica fields without hiding boundary " + "or applicability uncertainty." + ), + "archetype": "single-panel quantitative classification figure", + "backend": "Python/matplotlib", + "role": "classification evidence", + "evidence": "SiO2 and Na2O + K2O in wt%", + "target_output": "publication figure", + "review_risks": [ + "non-volcanic samples", + "volatile-free normalization not confirmed", + "alteration and alkali mobility", + "boundary-point ambiguity", + "unresolved normative subtypes", + ], + }, + "configuration": { + "sample_column": sample_column, + "group_column": group_column, + "groups": groups, + "confirm_volcanic": confirm_volcanic, + "composition_basis": composition_basis, + "x": "SiO2", + "y": "Na2O + K2O", + "unit": "wt%", + "fixed_model_limits": True, + "legend_layout": legend_layout, + "legend_position": plot_info["legend_position"], + "width_mm": width_mm, + "height_mm": height_mm, + "png_dpi": dpi, + "tiff_dpi": dpi, + }, + "plot": plot_info, + "outputs": [ + shareable_file_record(path) for path in figure_paths + ], + "source_data": shareable_file_record(source_path), + "submission_qa": { + "final_size_mm": [width_mm, height_mm], + "svg_text_editable": True, + "pdf_font_type": 42, + "raster_dpi": dpi, + "tiff_compression": "LZW", + "white_background": True, + "classification_boundaries_versioned": True, + "colourblind_support": "group colour plus marker shape", + "source_data_exported": True, + }, + "report_file": report_path.name, + "issues": run_issues, + "interpretation_guidance": [ + "Report the declared composition basis and the TAS reference in the caption or methods.", + "Treat boundary and outside-model analyses as requiring manual review.", + "Do not use TAS alone to establish petrogenesis, magma series, or tectonic setting.", + "Review alteration and alkali mobility before accepting a volcanic rock name.", + "Resolve Trachyte/Trachydacite and Tephrite/Basanite only with the additional normative criteria.", + ], + "scientific_caveat": ( + "TAS is a chemical nomenclature tool for its declared volcanic " + "domain; it is not a stand-alone genetic or tectonic classifier." + ), + } + report_path.write_text( + json.dumps(report, ensure_ascii=False, indent=2) + "\n", + encoding="utf-8", + ) + return report + except ( + InspectionError, + PlottingError, + OSError, + ValueError, + ) as exc: + return tas_error(input_path, str(exc)) + finally: + if figure is not None: + plt.close(figure) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="生成带 IUGS 火山岩分类字段的投稿级 TAS 图。" + ) + parser.add_argument("input", type=Path, help="CSV、TXT 或 Excel 输入表格") + parser.add_argument("--output-dir", type=Path, required=True) + parser.add_argument("--stem", help="不含扩展名的输出文件名") + parser.add_argument("--sheet", help="Excel 工作表名称或编号") + parser.add_argument("--sample-column", help="明确指定样品编号列") + parser.add_argument("--group-column", help="明确指定可选分组列") + parser.add_argument("--groups", help="逗号分隔的待绘制分组") + parser.add_argument( + "--confirm-volcanic", + action="store_true", + help="明确确认样品属于火山岩 TAS 适用范围", + ) + parser.add_argument( + "--composition-basis", + choices=COMPOSITION_BASES[1:], + help="声明输入为 anhydrous-normalized 或 as-reported", + ) + parser.add_argument("--title", help="可选图题") + parser.add_argument("--width-mm", type=float, default=150.0) + parser.add_argument("--height-mm", type=float, default=120.0) + parser.add_argument("--dpi", type=int, default=600) + parser.add_argument( + "--legend-layout", + choices=LEGEND_LAYOUTS, + default="inside-auto", + ) + parser.add_argument("--overwrite", action="store_true") + return parser.parse_args() + + +def main() -> int: + if hasattr(sys.stdout, "reconfigure"): + sys.stdout.reconfigure(encoding="utf-8") + if hasattr(sys.stderr, "reconfigure"): + sys.stderr.reconfigure(encoding="utf-8") + args = parse_args() + report = plot_tas_path( + args.input, + args.output_dir, + stem=args.stem, + requested_sheet=args.sheet, + requested_sample_column=args.sample_column, + requested_group_column=args.group_column, + requested_groups=args.groups, + confirm_volcanic=args.confirm_volcanic, + composition_basis=args.composition_basis or "unknown", + title=args.title, + width_mm=args.width_mm, + height_mm=args.height_mm, + dpi=args.dpi, + legend_layout=args.legend_layout, + overwrite=args.overwrite, + ) + json.dump(report, sys.stdout, ensure_ascii=False, indent=2) + sys.stdout.write("\n") + print(f"TAS 绘图完成：{report['status']}", file=sys.stderr) + if report["status"] == "ready": + return 0 + if report["status"] == "error": + return 1 + return 2 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/test_inspect_major_data.py b/tests/test_inspect_major_data.py new file mode 100644 index 0000000..27ad927 --- /dev/null +++ b/tests/test_inspect_major_data.py @@ -0,0 +1,146 @@ +import json +import subprocess +import sys +from pathlib import Path + +import numpy as np +import pandas as pd + + +ROOT = Path(__file__).resolve().parents[1] +SKILL = ROOT / "skills" / "geoskills" +SCRIPT = SKILL / "scripts" / "inspect_major_data.py" +EXAMPLE = SKILL / "examples" / "synthetic_major_element_data.csv" +SCRIPTS = SKILL / "scripts" +sys.path.insert(0, str(SCRIPTS)) + +from inspect_major_data import ( # noqa: E402 + inspect_major_frame, + prepare_geochem_frame, + read_major_table, +) + + +def run_inspector( + *arguments: object, +) -> tuple[subprocess.CompletedProcess[str], dict]: + result = subprocess.run( + [sys.executable, str(SCRIPT), *(str(item) for item in arguments)], + check=False, + capture_output=True, + text=True, + encoding="utf-8", + ) + return result, json.loads(result.stdout) + + +def test_valid_flat_major_csv_is_ready_and_private() -> None: + result, report = run_inspector(EXAMPLE) + + assert result.returncode == 0 + assert report["status"] == "ready" + assert report["shape"] == {"rows": 10, "columns": 15} + assert report["sample_id_candidates"] == ["Sample"] + assert report["group_candidates"] == ["Group"] + recognized = { + item["analyte"]: item + for item in report["analytes"]["recognized"] + } + assert recognized["SiO2"]["unit"] == "wt%" + assert recognized["Rb"]["unit"] == "ppm" + assert "filename" not in report["source"] + assert "path" not in report["source"] + + +def test_transposed_table_tracks_major_and_trace_sections( + tmp_path: Path, +) -> None: + path = tmp_path / "supplement.xlsx" + pd.DataFrame( + [ + ["Rock type", "Suite A", None], + ["Sample No.", "S1", "S2"], + ["Major element (wt.%)", None, None], + ["SiO2", 50.0, 52.0], + ["MgO", 6.0, 4.5], + ["Na2O", 3.0, 3.4], + ["K2O", 1.0, 1.4], + ["Trace element (ppm)", None, None], + ["Rb", 35, 48], + ["Zr", 120, 145], + ] + ).to_excel(path, index=False, header=False) + + frame, source = read_major_table(path, None) + assert frame is not None + report = inspect_major_frame(frame, source) + + assert report["status"] == "ready" + assert source["layout"] == "column_per_sample_transposed" + assert list(frame["Group"]) == ["Suite A", "Suite A"] + units = { + item["analyte"]: item["unit"] + for item in report["analytes"]["recognized"] + } + assert units["SiO2"] == "wt%" + assert units["Rb"] == "ppm" + + +def test_unknown_major_unit_requires_review() -> None: + frame = pd.DataFrame( + { + "Sample": ["S1", "S2"], + "SiO2": [50, 52], + "MgO": [5, 4], + } + ) + report = inspect_major_frame(frame, {"format": ".csv"}) + + assert report["status"] == "needs_review" + assert sum(item["code"] == "E621" for item in report["issues"]) == 2 + + +def test_negative_concentration_requires_review() -> None: + frame = pd.DataFrame( + { + "Sample": ["S1", "S2"], + "SiO2_wt%": [50, 52], + "MgO_wt%": [5, -0.1], + } + ) + report = inspect_major_frame(frame, {"format": ".csv"}) + + assert report["status"] == "needs_review" + assert any(item["code"] == "E632" for item in report["issues"]) + + +def test_bdl_trace_value_is_preserved_as_missing() -> None: + frame = pd.DataFrame( + { + "Sample": ["S1", "S2"], + "SiO2_wt%": [50, 52], + "MgO_wt%": [5, 4], + "Rb_ppm": ["bdl", 20], + } + ) + report = inspect_major_frame(frame, {"format": ".csv"}) + canonical = prepare_geochem_frame(frame, report) + + assert report["status"] == "ready" + assert any(item["code"] == "W631" for item in report["issues"]) + assert np.isnan(canonical.at[0, "Rb"]) + + +def test_duplicate_canonical_analyte_requires_review() -> None: + frame = pd.DataFrame( + { + "Sample": ["S1", "S2"], + "SiO2_wt%": [50, 52], + "silica_wt%": [50, 52], + "MgO_wt%": [5, 4], + } + ) + report = inspect_major_frame(frame, {"format": ".csv"}) + + assert report["status"] == "needs_review" + assert any(item["code"] == "E611" for item in report["issues"]) diff --git a/tests/test_plot_harker.py b/tests/test_plot_harker.py new file mode 100644 index 0000000..d40adb1 --- /dev/null +++ b/tests/test_plot_harker.py @@ -0,0 +1,211 @@ +import json +import subprocess +import sys +from pathlib import Path + +import pandas as pd +import pytest +from PIL import Image + + +ROOT = Path(__file__).resolve().parents[1] +SKILL = ROOT / "skills" / "geoskills" +SCRIPT = SKILL / "scripts" / "plot_harker.py" +EXAMPLE = SKILL / "examples" / "synthetic_major_element_data.csv" + + +def run_plotter( + *arguments: object, +) -> tuple[subprocess.CompletedProcess[str], dict]: + result = subprocess.run( + [sys.executable, str(SCRIPT), *(str(item) for item in arguments)], + check=False, + capture_output=True, + text=True, + encoding="utf-8", + ) + return result, json.loads(result.stdout) + + +def test_exports_harker_publication_bundle(tmp_path: Path) -> None: + output = tmp_path / "harker" + result, report = run_plotter( + EXAMPLE, + "--output-dir", + output, + "--stem", + "harker_test", + "--y", + "MgO,CaO,K2O", + "--width-mm", + "100", + "--height-mm", + "70", + "--dpi", + "100", + ) + paths = { + suffix: output / f"harker_test.{suffix}" + for suffix in ("svg", "pdf", "tiff", "png") + } + + assert result.returncode == 0 + assert report["status"] == "ready" + assert report["configuration"]["x"] == "SiO2" + assert report["configuration"]["y"] == ["MgO", "CaO", "K2O"] + assert report["plot"]["legend_position"] == "shared_figure_top" + assert report["plot"]["legend_within_figure"] is True + assert report["plot"]["shared_x_label"] is True + assert report["plot"]["axes_frame"] == "full" + assert all(path.exists() and path.stat().st_size > 0 for path in paths.values()) + assert " None: + output = tmp_path / "custom" + result, report = run_plotter( + EXAMPLE, + "--output-dir", + output, + "--stem", + "custom", + "--x", + "MgO", + "--y", + "Rb,Zr", + "--groups", + "Series A,Series B", + "--dpi", + "90", + ) + + assert result.returncode == 0 + assert report["configuration"]["x"] == "MgO" + assert report["configuration"]["y"] == ["Rb", "Zr"] + assert report["configuration"]["groups"] == ["Series A", "Series B"] + assert report["plot"]["sample_count"] == 6 + + +def test_harker_limits_are_clean_and_do_not_clip() -> None: + sys.path.insert(0, str(SKILL / "scripts")) + from plot_geochem_common import clean_linear_limits + from plot_harker import ( + automatic_columns, + potentially_mixed_harker_groups, + ) + + limits = clean_linear_limits(pd.Series([7.89, 1437.96])) + + assert limits["lower"] <= 7.89 + assert limits["upper"] >= 1437.96 + assert float(limits["lower"]).is_integer() + assert float(limits["upper"]).is_integer() + assert automatic_columns(8) == 4 + assert automatic_columns(6) == 3 + assert automatic_columns(4) == 2 + assert automatic_columns(3) == 3 + assert automatic_columns(9) == 3 + assert potentially_mixed_harker_groups( + ["Syenite", "Carbonatite suite"] + ) == ["Carbonatite suite"] + assert potentially_mixed_harker_groups(["Carbonatite A"]) == [] + + +def test_missing_pair_is_reported_not_imputed(tmp_path: Path) -> None: + input_path = tmp_path / "missing.csv" + output = tmp_path / "missing" + frame = pd.read_csv(EXAMPLE) + frame.loc[0, "MgO_wt%"] = None + frame.to_csv(input_path, index=False) + + result, report = run_plotter( + input_path, + "--output-dir", + output, + "--stem", + "missing", + "--y", + "MgO,CaO", + "--dpi", + "90", + ) + + assert result.returncode == 0 + assert report["plot"]["missing_pairs"]["MgO"] == 1 + assert any(item["code"] == "W651" for item in report["issues"]) + + +def test_existing_bundle_requires_explicit_overwrite( + tmp_path: Path, +) -> None: + output = tmp_path / "existing" + first, _ = run_plotter( + EXAMPLE, + "--output-dir", + output, + "--stem", + "same", + "--y", + "MgO,CaO", + "--dpi", + "90", + ) + second, report = run_plotter( + EXAMPLE, + "--output-dir", + output, + "--stem", + "same", + "--y", + "MgO,CaO", + "--dpi", + "90", + ) + + assert first.returncode == 0 + assert second.returncode == 1 + assert report["status"] == "error" + + +def test_long_group_legend_reflows_inside_figure(tmp_path: Path) -> None: + input_path = tmp_path / "long_groups.csv" + output = tmp_path / "long_groups" + frame = pd.read_csv(EXAMPLE) + frame["Group"] = [ + "Long alkali feldspar suite", + "Long alkali feldspar suite", + "Long nepheline-bearing suite", + "Long nepheline-bearing suite", + "Long evolved syenitic suite", + "Long evolved syenitic suite", + "Long mixed-lithology review suite", + "Long mixed-lithology review suite", + "Long mixed-lithology review suite", + "Long mixed-lithology review suite", + ] + frame.to_csv(input_path, index=False) + + result, report = run_plotter( + input_path, + "--output-dir", + output, + "--stem", + "long_groups", + "--y", + "MgO,CaO,K2O,P2O5", + "--dpi", + "90", + ) + + assert result.returncode == 0 + assert report["plot"]["legend_columns"] == 2 + assert report["plot"]["legend_rows"] == 2 + assert report["plot"]["legend_within_figure"] is True diff --git a/tests/test_plot_tas.py b/tests/test_plot_tas.py new file mode 100644 index 0000000..5863cfc --- /dev/null +++ b/tests/test_plot_tas.py @@ -0,0 +1,178 @@ +import json +import subprocess +import sys +from pathlib import Path + +import pytest +from PIL import Image + + +ROOT = Path(__file__).resolve().parents[1] +SKILL = ROOT / "skills" / "geoskills" +SCRIPT = SKILL / "scripts" / "plot_tas.py" +EXAMPLE = SKILL / "examples" / "synthetic_major_element_data.csv" +MODEL = SKILL / "assets" / "classification" / "tas-lemaitre-2002.json" +sys.path.insert(0, str(SKILL / "scripts")) + +from plot_tas import classify_tas_point, load_tas_model # noqa: E402 + + +def run_plotter( + *arguments: object, +) -> tuple[subprocess.CompletedProcess[str], dict]: + result = subprocess.run( + [sys.executable, str(SCRIPT), *(str(item) for item in arguments)], + check=False, + capture_output=True, + text=True, + encoding="utf-8", + ) + return result, json.loads(result.stdout) + + +def test_tas_asset_has_expected_fields_and_vertices() -> None: + model = load_tas_model(MODEL) + fields = {field["id"]: field for field in model["fields"]} + + assert model["source"]["boundary_doi"] == "10.1007/BF01160698" + assert len(fields) == 15 + assert fields["B"]["polygon"] == [ + [45.0, 0.0], + [45.0, 5.0], + [52.0, 5.0], + [52.0, 0.0], + ] + assert fields["T"]["name"] == "Trachyte/Trachydacite" + assert fields["U1"]["name"] == "Tephrite/Basanite" + + +@pytest.mark.parametrize( + ("x", "y", "field"), + [ + (47.0, 3.0, "B"), + (54.0, 4.0, "O1"), + (59.5, 4.8, "O2"), + (66.0, 5.5, "O3"), + (74.0, 6.2, "R"), + (50.0, 6.0, "S1"), + (54.0, 7.2, "S2"), + (59.0, 9.0, "S3"), + (49.0, 10.3, "U2"), + (58.0, 14.0, "Ph"), + ], +) +def test_representative_points_classify_deterministically( + x: float, + y: float, + field: str, +) -> None: + result = classify_tas_point(x, y, load_tas_model(MODEL)) + + assert result["status"] == "classified" + assert result["field"] == field + + +def test_point_on_shared_boundary_requires_review() -> None: + result = classify_tas_point(52.0, 5.0, load_tas_model(MODEL)) + + assert result["status"] == "boundary_review" + assert result["name"] is None + + +def test_tas_requires_domain_and_basis_confirmation( + tmp_path: Path, +) -> None: + result, report = run_plotter( + EXAMPLE, + "--output-dir", + tmp_path / "blocked", + ) + + assert result.returncode == 2 + assert report["status"] == "blocked" + assert {item["code"] for item in report["issues"]} == { + "E711", + "E712", + } + + +def test_exports_tas_classification_bundle(tmp_path: Path) -> None: + output = tmp_path / "tas" + result, report = run_plotter( + EXAMPLE, + "--output-dir", + output, + "--stem", + "tas_test", + "--confirm-volcanic", + "--composition-basis", + "anhydrous-normalized", + "--width-mm", + "100", + "--height-mm", + "80", + "--dpi", + "100", + ) + paths = { + suffix: output / f"tas_test.{suffix}" + for suffix in ("svg", "pdf", "tiff", "png") + } + + assert result.returncode == 0 + assert report["status"] == "ready" + assert report["model"]["id"] == "TAS_LeMaitre2002_Volcanic_CombinedT" + assert report["plot"]["classification_status_counts"] == { + "classified": 10 + } + assert report["plot"]["axes_frame"] == "full" + assert report["configuration"]["legend_position"].startswith("inside_") + assert "filename" not in report["source"] + assert "path" not in report["source"] + svg = paths["svg"].read_text(encoding="utf-8") + assert " None: + result, report = run_plotter( + EXAMPLE, + "--output-dir", + tmp_path / "reported", + "--confirm-volcanic", + "--composition-basis", + "as-reported", + "--dpi", + "90", + ) + + assert result.returncode == 0 + assert any(item["code"] == "W711" for item in report["issues"]) + + +def test_existing_tas_bundle_requires_overwrite(tmp_path: Path) -> None: + output = tmp_path / "existing" + arguments = ( + EXAMPLE, + "--output-dir", + output, + "--confirm-volcanic", + "--composition-basis", + "anhydrous-normalized", + "--dpi", + "90", + ) + first, _ = run_plotter(*arguments) + second, report = run_plotter(*arguments) + + assert first.returncode == 0 + assert second.returncode == 1 + assert report["status"] == "error" From 67ec8f84f6bc9cbda283caa2bef1ed85eb8da851 Mon Sep 17 00:00:00 2001 From: gronbow <59046437+gronbow@users.noreply.github.com> Date: Thu, 30 Jul 2026 18:36:57 +0800 Subject: [PATCH 2/2] Prepare GeoSkills v0.3.0 release --- AGENTS.md | 4 ++-- CHANGELOG.md | 20 ++++++++++++++++++++ README.md | 20 ++++++++++++-------- skills/geoskills/SKILL.md | 8 ++++---- skills/geoskills/references/data-contract.md | 6 +++--- 5 files changed, 41 insertions(+), 17 deletions(-) create mode 100644 CHANGELOG.md diff --git a/AGENTS.md b/AGENTS.md index 146977c..e7096e4 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -9,6 +9,6 @@ - Write machine-readable results to standard output and concise diagnostics to standard error where practical. - Add tests for numerical rules, invalid input, missing values, non-positive values, and exported files. - Do not claim the skill is research-ready until scientific fixtures and end-to-end tests pass. -- Preserve the published v0.1.0 REE workflow and the reviewed v0.2.0 spider workflow. -- GeoSkills v0.3.0 may add Harker variation diagrams and volcanic TAS classification only; defer isotope, tectonic-discrimination, and other diagram families to later reviewed versions. +- Preserve the published v0.3.0 REE, spider, Harker, and volcanic TAS workflows. +- Add post-v0.3.0 diagram families only through separate scientific review and validation; defer isotope and tectonic-discrimination diagrams until those workflows are explicitly reviewed. - Require explicit volcanic applicability and composition-basis declarations before TAS plotting; never silently classify intrusive or otherwise out-of-scope samples. diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..ebf999e --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,20 @@ +# Changelog + +This file records the user-visible changes in each public GeoSkills release. + +## v0.3.0 — 2026-07-30 + +The reviewed v0.2.0 spider-diagram milestone was not tagged separately; it is first published as part of this complete v0.3.0 release. + +- Added validated primitive-mantle and N-MORB normalized trace-element spider diagrams. +- Added customizable multi-panel Harker variation diagrams with clean outward-rounded axes and shared group legends. +- Added guarded volcanic TAS classification with versioned Le Maitre/Le Bas boundaries, explicit composition-basis declarations, and boundary-review states. +- Added flat and unambiguous transposed CSV, TXT, and XLSX inspection for the new workflows. +- Added editable SVG/PDF, 600 dpi PNG/TIFF, plotted-data CSV, and machine-readable JSON export bundles. +- Preserved below-detection-limit values as missing and kept private inputs and generated outputs outside Git. +- Validated the complete release with 69 automated tests and Ubuntu/Windows CI on Python 3.11 and 3.12. + +## v0.1.0 — 2026-07-28 + +- Published the first GeoSkills release for Sun and McDonough (1989) C1 chondrite-normalized REE patterns. +- Added local table inspection, deterministic normalization, submission-oriented exports, and scientific/privacy safeguards. diff --git a/README.md b/README.md index 1ea1570..81c3758 100644 --- a/README.md +++ b/README.md @@ -2,8 +2,7 @@ GeoSkills 是一个面向地质学与地球化学研究的本地 Agent Skill。 -公开的 v0.1.0 只提供球粒陨石标准化 REE 配分图。 -v0.2.0 蛛网图工作流已经完成人工审核。当前本地 `feature/harker-tas-v0.3` 分支是 v0.3.0 审核候选版：保留前两类图件，并新增 Harker 变化图和带适用性保护的火山岩 TAS 分类图。该分支在用户审核前不会上传 GitHub。 +当前正式版本为 [v0.3.0](https://github.com/gronbow/GeoSkills/releases/tag/v0.3.0)，包含球粒陨石标准化 REE 配分图、微量元素蛛网图、Harker 变化图，以及带适用性保护的火山岩 TAS 分类图。 ## 当前可以做什么 @@ -23,27 +22,31 @@ v0.2.0 蛛网图工作流已经完成人工审核。当前本地 `feature/harker ## 安装到 Codex -以下公开仓库的默认分支目前仍是经过审核的 v0.1.0 REE 版本；v0.2.0 蛛网图通过功能分支和 Pull Request 接受合并检查。v0.3.0 的 Harker/TAS 代码仍处于本地审核阶段： +公开仓库的默认分支 `main` 包含经过审核和自动测试的正式 v0.3.0： ```text https://github.com/gronbow/GeoSkills ``` -也可以手动把 `skills/geoskills` 文件夹复制到个人 Skill 目录： +也可以手动把 `skills/geoskills` 文件夹复制到 Codex 的个人 Skill 目录： ```text -%USERPROFILE%\.agents\skills\geoskills +%USERPROFILE%\.codex\skills\geoskills ``` +如果设置了 `CODEX_HOME`，则复制到 `%CODEX_HOME%\skills\geoskills`。复制完成后，重新开启一个 Codex 任务即可使用。 + ## Windows 快速开始 -在 PowerShell 中进入项目目录，安装依赖： +需要 Python 3.11 或 3.12。先运行 `python --version` 确认版本，然后在 PowerShell 中进入项目目录并安装依赖： ```powershell -py -m venv .venv +python -m venv .venv .\.venv\Scripts\python.exe -m pip install -r requirements-dev.txt ``` +如果系统只提供 Python Launcher，可以把第一条命令改为 `py -3.12 -m venv .venv`。 + 先检查示例数据： ```powershell @@ -106,6 +109,7 @@ TAS 命令要求明确确认样品属于火山岩，并声明数据是无水归 ```text GeoSkills/ ├── README.md +├── CHANGELOG.md ├── AGENTS.md ├── requirements-dev.txt ├── skills/ @@ -135,7 +139,7 @@ GeoSkills/ ## 当前状态 -GeoSkills v0.1.0 是仅包含 REE 配分图的公开测试版。v0.2.0 蛛网图发布候选版已通过 41 项自动测试、真实数据测试、导出审计和人工图形复核。v0.3.0 Harker/TAS 候选版已在本地通过 69 项自动测试、已发表数据的 Harker 验证、合成火山岩 TAS 验证、导出审计和人工图面复核；用户确认前不会提交或上传。 +GeoSkills v0.3.0 是当前正式版本。它整合了 v0.1.0 的 REE 工作流、经过审核的微量元素蛛网图，以及新增的 Harker 和火山岩 TAS 工作流。完整版本通过了 69 项自动测试、已发表数据的 Harker 验证、合成火山岩 TAS 验证、导出审计、人工图面复核，以及 Ubuntu/Windows 上 Python 3.11/3.12 的 GitHub Actions 检查。版本变化见 [CHANGELOG.md](CHANGELOG.md)。 ## 许可 diff --git a/skills/geoskills/SKILL.md b/skills/geoskills/SKILL.md index cb08433..136a9be 100644 --- a/skills/geoskills/SKILL.md +++ b/skills/geoskills/SKILL.md @@ -3,13 +3,13 @@ name: geoskills description: "Inspect whole-rock geochemical tables and create validated, submission-oriented REE patterns, trace-element spider diagrams, Harker variation diagrams, or volcanic TAS classification plots from CSV, TXT, or Excel data. Use GeoSkills to identify sample, group, element, and oxide columns; verify ppm or wt% units and invalid values; apply versioned Sun and McDonough (1989) references; require explicit TAS applicability declarations; export editable SVG/PDF plus 600 dpi TIFF/PNG and exact plotted data; or explain why a table cannot yet be plotted safely. Do not use the current version for isotope, tectonic-discrimination, or other diagram families." --- -# GeoSkills v0.3.0 review candidate +# GeoSkills v0.3.0 Create reproducible REE patterns, trace-element spider diagrams, Harker variation diagrams, and volcanic TAS plots through deterministic local Python scripts. Use the language model to select and explain the workflow, never to calculate normalized ratios or classification fields manually. -## Development status +## Release status -Treat v0.1.0 REE plotting as the stable published baseline. The v0.2.0 spider workflow has passed scientific fixtures, real-data output, export QA, Skill validation, and user review. Harker and TAS are v0.3.0 review-candidate workflows and must not be described as published until their local review and release process are complete. Do not imply support for isotope or tectonic-discrimination diagrams. +Treat v0.3.0 as the current released workflow for REE patterns, trace-element spider diagrams, Harker variation diagrams, and guarded volcanic TAS classification. These workflows have passed scientific fixtures, real-data or synthetic-data validation as appropriate, export QA, Skill validation, cross-platform automated tests, and user review. Do not imply support for isotope or tectonic-discrimination diagrams. ## Route the request @@ -70,7 +70,7 @@ TAS uses `SiO2` versus `Na2O + K2O` in wt%. Never infer volcanic applicability f - Convert only explicit `K2O`, `P2O5`, and `TiO2` wt% columns to K, P, and Ti ppm. Record the CIAAW/IUPAC atomic weights, formula, and factor used. - Reject finite zero and negative values on logarithmic axes. - Preserve the cited Sun and McDonough (1989) incompatibility order even when the user selects a subset. -- Use only the versioned local assets `PrimitiveMantle_SM89`, `PrimitiveMantleModified_SM89`, and `NMORB_SM89` for the spider draft. +- Use only the versioned local assets `PrimitiveMantle_SM89`, `PrimitiveMantleModified_SM89`, and `NMORB_SM89` for the spider workflow. - Use the source footnote's modified primitive mantle as the spider default. Keep the printed and modified variants separate; do not silently replace Cs or Pb, and warn when an affected element is plotted with the printed variant. - Set log limits from finite positive ratios, add a declared margin, round to clean decimal bounds without clipping, and show unity only when it lies inside the range. - Use the versioned `TAS_LeMaitre2002_Volcanic_CombinedT` asset; keep `Trachyte/Trachydacite` and `Tephrite/Basanite` unresolved when normative information is absent. diff --git a/skills/geoskills/references/data-contract.md b/skills/geoskills/references/data-contract.md index b76d85f..e615d1b 100644 --- a/skills/geoskills/references/data-contract.md +++ b/skills/geoskills/references/data-contract.md @@ -1,12 +1,12 @@ # REE input data contract -This file records the first-release rules for data accepted by `geoskills` v1. +This file records the REE input rules introduced in GeoSkills v0.1.0 and retained in v0.3.0. ## Required information - One sample identifier column with a non-empty value for every row. - REE concentration columns selected from `La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, Lu`. -- A confirmed concentration unit. The first release accepts ppm. +- A confirmed concentration unit. Direct REE concentrations are accepted in ppm. ## Optional information @@ -30,7 +30,7 @@ This file records the first-release rules for data accepted by `geoskills` v1. - Preserve blank cells as missing values. - Treat strings such as `<0.01`, `BDL`, and `n.d.` as detection-limit states, not as zero. - Reject zero and negative values for logarithmic REE plots. -- Do not interpolate missing REE concentrations in the first release. +- Do not interpolate missing REE concentrations. ## File handling