# geovector-3d-mcp

**Downhole 3-D geochemical vectoring with honest uncertainty, plus an MCP server so AI
agents can call the maths instead of improvising it.**

You have six holes through a pegmatite, all sub-economic. The fractionation indices
change from hole to hole: K/Rb falls, Nb/Ta falls, Cs rises. The question every
exploration team asks next is which way, and how sure are we? This package answers
it from collar, survey, assay and lithology tables:

1. **Desurvey** every hole by minimum curvature, so sample positions are true 3-D
   points, not collar-plus-depth guesses.
2. **Clean the geochemistry.** It parses below-detection values, uses log-ratios
   instead of raw ratios, and computes indices only in the rock they are defined for.
3. **Fit a local 3-D gradient** of the index, aligned to the host structure, with a
   95 % confidence cone on its direction and a significance test that it exists.
4. **Propose a step-out hole** (collar, azimuth, dip, target depth), or decline to
   when the direction is not resolved.

![Synthetic benchmark](docs/benchmark_3d.png)

> The benchmark data in this repository are **synthetic**, generated by
> [`scripts/make_benchmark.py`](scripts/make_benchmark.py) with a planted answer.
> Recovering that answer shows the code is correct. It does not show that real
> pegmatites zone this simply.

---

## Quick start

```bash
git clone https://github.com/YASHUJOSHI/geovector-3d-mcp.git
cd geovector-3d-mcp
pip install -e ".[dev]"          # core + MCP + Streamlit + pyproj + pytest
pytest                           # 43 tests
streamlit run app/streamlit_viewer.py
```

From Python:

```python
from geovector3d import run_vectoring

r = run_vectoring("data/benchmark_stuck_block", index="ln_K_Rb",
                  host_codes=("PEG",), strike=300, dip=50)
r["gradient"]["direction_trend_deg"], r["gradient"]["cone_half_angle_deg"]
r["proposal"]          # collar, azimuth, dip, target depth, confidence, caveats
```

### Input tables

| File | Required columns | Notes |
|---|---|---|
| `collar.csv` | `BHID, X, Y, Z, TOTAL_DEPTH` | local grid or projected CRS, metres |
| `survey.csv` | `BHID, DEPTH, AZIMUTH, DIP` | dip **negative down**; missing 0 m station is copied from the first shot |
| `assay.csv` | `BHID, FROM, TO, <analytes…>` | `<0.5` or `-0.5` both mean "below 0.5" |
| `lithology.csv` | `BHID, FROM, TO, ROCK_CODE` | optional but strongly recommended |
| `units.json` | `{analyte: "ppm" \| "ppb" \| "pct" \| "oxide_pct"}` | optional; default is ppm |

## Results on the benchmark

Six inclined holes cut a 20 m pegmatite dyke (strike 300°, dip 50° NE). K/Rb falls
from 140 to 38 and Nb/Ta from 12 to 2.5 down a 34° plunge toward 335°. Planted truth:
trend 335.0° / plunge 34.4°.

| Configuration | Recovered trend / plunge | Error | 95 % cone | Outcome |
|---|---|---|---|---|
| ln K/Rb, pegmatite only, structure-aligned | 333.6° / 33.4° | 1.5° | 10.5° | step-out proposed, high confidence |
| ln Nb/Ta, same | | 7.5° | 10.9° | proposed |
| ln Zr/Hf, same (weaker planted trend) | | 11.6° | 37.3° | moderate |
| ln K/Rb, isotropic search (no structure) | | 38° | 45° | borderline |
| ln K/Rb, **no lithology filter** | | | | gradient not significant (Wald p = 0.43): no recommendation |

The last two rows are the point of the package. Without structure, and without
keeping wall rock out of the fit, the same data give a wrong or empty answer. A tool
that reported a direction regardless would have sent a hole the wrong way.

---

## Method

[`METHODOLOGY.md`](METHODOLOGY.md) has the derivations. In short:

### Minimum-curvature desurveying

Each survey station becomes a unit tangent
**t** = (cos δ sin α, cos δ cos α, sin δ), with azimuth α and dip δ negative down.
Between stations *i* and *i+1*, separated by ΔMD:

- dogleg angle β = 2 arcsin(|**t**ᵢ₊₁ − **t**ᵢ| / 2)
- ratio factor RF = (2/β) tan(β/2)
- displacement Δ**x** = (ΔMD/2) · RF · (**t**ᵢ + **t**ᵢ₊₁)

β is computed from the chord rather than arccos(**t**ᵢ·**t**ᵢ₊₁), which loses most of
its precision for the nearly straight courses that dominate real surveys. RF is
replaced by its limit 1 + β²/12 when β < 10⁻⁶. A sample midpoint between stations is
placed on the arc, not the chord. The sub-arc to a fraction *f* ends at the slerp of
the two tangents and is itself a minimum-curvature course of length *f*·ΔMD. The
tests check this against a constant-build circular arc to 10⁻⁸ m.

Because everything is carried as vectors, a vertical hole is not a special case. At
δ = −90° the azimuth multiplies zero and drops out, so there is no gimbal lock.

### Compositional data

- **Non-detects** are kept as censored values and replaced multiplicatively
  (Martín-Fernández et al., 2003) at 0.65 × LLD. This leaves ratios among detected
  parts unchanged. Analytes censored in more than 30 % of samples are reported as
  unusable, and a ratio of two censored values is set to missing.
- **Vectoring runs on log-ratios** (ln K/Rb, ln Nb/Ta, …). A pairwise log-ratio is
  scale-invariant and subcompositionally coherent. Raw ratios are skewed, and their
  gradients are dominated by the high end. CLR and ILR are provided for
  multi-element work but are computed on a fixed, closed composition, because CLR
  scores change whenever the analyte suite changes.
- **Indices are host-restricted.** LCT ratios are computed only in pegmatite and
  granite codes, and the Ishikawa AI and CCPI only in felsic volcanics.

### Local gradient

A weighted least-squares fit of f(**x**) = a + **g**·(**x** − **x**₀) around **x**₀,
with three safeguards for the way drillholes sample space:

| Problem | Safeguard |
|---|---|
| Dense down-hole, sparse between holes, so the fit is driven by along-hole variation | at most 3 samples per hole per fit |
| Isotropic search crosses barren contacts | distances measured in a strike / dip / pole frame with separate ranges; optionally fit only the in-plane components |
| A direction is always returned, even from noise | a Wald test that **g** ≠ 0, and a 95 % cone half-angle from the parameter covariance |

The cone is arctan(√(χ²₂,₀.₉₅ · λ⊥) / |**g**|), where λ⊥ is the largest variance
perpendicular to **g**. It is a better guide than the design-matrix condition
number, which is reported too, because it scales with the signal.

### Step-out design

The target is placed along the in-plane fractionation direction, starting from the
most-evolved decile of host samples. The default step is the data-support distance
(median spacing of host intercepts). The hole is aimed down the host pole to cut the
body at a high angle, with dip clipped to [−85°, −45°]. Longer steps are allowed but
flagged `beyond_support`. A local gradient cannot tell you whether a trend continues
for another 120 m: that is the hypothesis the hole tests.

### Is a barren surface sample barren?

`diagnose_false_surface_anomaly` checks that immobile ratios (Zr/Hf, Nb/Ta, Ti/Zr)
match the fresh parent, which rules out transported cover. It then computes tau
mass-transfer coefficients (Brimhall & Dietrich, 1987) for mobile elements. It
returns evidence and caveats, not a verdict, because Zr and Ti are only relatively
immobile.

---

## MCP server for AI agents

Agents should not do drillhole trigonometry or log-ratio algebra in their heads. The
server exposes the engine as four deterministic tools. Every response states its sign
conventions, and bad input returns an explicit error.

| Tool | Does |
|---|---|
| `validate_and_desurvey_drillholes` | checks tables, desurveys, reports end-of-hole positions, dogleg severity and copied collar stations |
| `compute_3d_fractionation_vector` | host-filtered log-ratio gradient with cone and significance |
| `propose_stepout_borehole` | collar / azimuth / dip / depth, or `no_recommendation` |
| `diagnose_false_surface_anomaly` | immobile-ratio check + tau mass transfer |

Built on the MCP Python SDK 2.x (`MCPServer`, formerly `FastMCP`). Add it to Claude
Desktop's `claude_desktop_config.json`:

```json
{
  "mcpServers": {
    "geovector3d": { "command": "geovector3d-mcp" }
  }
}
```

If the console script is not on your PATH, use
`"command": "python", "args": ["-m", "geovector3d.mcp_server"]`.

## Repository layout

```
src/geovector3d/
  desurvey.py        minimum curvature, validation, interval mid-points
  compositional.py   censoring, units, CLR/ILR, LCT and alteration indices
  vector_field.py    structural frame, local gradient, cone, step-out design
  anomaly.py         immobile ratios and tau mass transfer
  workflow.py        load -> prepare -> gradient -> proposal
  export.py          CSV drill plan, GeoJSON (needs an EPSG code)
  mcp_server.py      MCP tools
app/streamlit_viewer.py
data/benchmark_stuck_block/   synthetic six-hole dataset + truth.json
scripts/                      benchmark and figure generators
tests/                        43 tests
```

## Limitations

- Linear local model: zoning that curves within the search radius is averaged.
- One structural orientation per run. Folded hosts need domaining; unfolding is out
  of scope.
- Fractionation indices assume a single magmatic system. Overprinting, alteration or
  multiple intrusive phases break the simple gradient picture.
- GeoJSON export needs a real CRS; local mine grids export as CSV.
- Validated only on the synthetic benchmark. Real-data case studies are the next step.

## How this was built

I wrote this code with an AI coding assistant. The geological framing is mine: the
choice of indices, host restriction, structure-aligned search, refusing to recommend
without a significant gradient, and flagging step-outs beyond data support.

## References

- Brimhall, G.H. & Dietrich, W.E. (1987). *Geochim. Cosmochim. Acta* 51, 567–587.
- Černý, P. (1991). Rare-element granitic pegmatites. *Geoscience Canada* 18, 49–67.
- Egozcue, J.J. et al. (2003). Isometric logratio transformations. *Math. Geol.* 35, 279–300.
- Ishikawa, Y. et al. (1976). *Mining Geology* 26, 105–117. Alteration index.
- Large, R.R. et al. (2001). *Econ. Geol.* 96, 957–971. AI and CCPI box plot.
- London, D. (2008). *Pegmatites.* Can. Mineral. Spec. Publ. 10.
- Martín-Fernández, J.A. et al. (2003). *Math. Geol.* 35, 253–278. Zero replacement.
- Sawaryn, S.J. & Thorogood, J.L. (2005). *SPE Drill. Complet.* 20, 77–86. Minimum curvature.

MIT licence.