geohashim

strata

A command-line tool that turns borehole logs into a layered model of the ground. For any point it gives the elevation of every contact with a 95% interval, and where that interval would be too wide to use, it refuses to answer.

A site investigation might drill forty holes over a square kilometre, so most of the ground has never been seen. Most tools still draw smooth, confident layers through it, and the section looks just as certain 150 m from a hole as it does next to one. strata says which parts are measured and which are guesses, and it checks its own intervals: of its 95% intervals, about 95% contain the truth.

View on GitHub How the intervals are checked

Plays the seven stages below in order, about 70 seconds. Scroll or press Esc to stop.

Language
Rust 1.74+, no dependencies
Tests
46 passing with cargo test, in a few seconds
Model
Ordinary kriging of the ground and each unit's thickness
Checked
95.3% and 94.0% inside the 95% intervals, 100 synthetic sites each
Refuses
Any contact with σ over 1.5 m, and exits 1
License
MIT
strata query examples/site.csv 490 440 at (490, 440), max sigma 1.5 m surface elevation 95% interval sigma ground surface 63.63 63.28 .. 63.97 0.17

All four are real runs of the release binary, output verbatim. examples/site.csv is the repo's 40-hole synthetic site, made by strata synth --seed 7. The first query is 4.9 m from hole BH16, the second 120.7 m from it, and BH16 is the nearest hole to both.

1Read the logs, and refuse to guess

The input is one CSV of interval logs: a units: line listing the column from top to bottom, then one row per logged interval, with depths below the collar. The parser fails closed. Anything ambiguous is an error that names the line, never a silent fix.

  • A misspelt unit, a gap or overlap between intervals, a unit logged twice or out of order, a hole that moves between rows, a hole's rows split up by another hole, and two holes at one spot.
  • A unit a hole skips between two units it did log is pinched out there: thickness 0, which is real information.
  • The unit a hole stops in has an unknown base. It is not used to fit anything, only checked afterwards by strata validate.

Every bad file exits 2. A model built on a typo would still look like a model, so a typo has to stop the run.

Eight bad files, eight line numbers

    
            
    
            
    typo.csv is the repo's own test file. The other seven are Illustrative input, written for this page in the same shape. Each error and exit code is what strata validate printed for that file.
    
          
    
        
    
      
    
    
    

    
    
      
    
        
    
          
    2Model thicknesses, so layers can't cross
    
          
    Interpolating each contact's elevation on its own lets layers cross: base of sand above base of clay, between two holes where neither does. strata interpolates the ground surface and each unit's thickness instead, clamps every thickness at zero, and stacks them.
    
          
      
            
    • Base of unit k is the ground minus every thickness down to k. Layers can't cross by construction, and a unit can thin out to nothing.
      • 
            
    • The variance of a contact is the sum of the variances above it, so uncertainty only grows downward.
      • 
            
    • The property tests look for a crossing at 400 random points on each of 25 random sites, smooth and rough, including up to 500 m outside the site. They find none.
      • 
          
    
          
    Adding the variances assumes the thickness errors are independent. Real thicknesses often aren't, which is one reason the intervals are measured in stage 4 instead of trusted.
    
        
    
        
    
          
    
            
    
              
    Two columns from strata query --json
    
              
            
    
            
    
              
    
              
    
            
    
            
    Elevations and σ are the query's own output at each point; each thickness is the difference between two of its elevations. At BH26 the log goes from topsoil straight to sand, so clay's thickness there is 0 and its σ is 0: the model passes through every hole exactly.
    
          
    
        
    
      
    
    
    

    
    
      
    
        
    
          
    3Krige each variable, with its error variance
    
          
    Each variable is interpolated by ordinary kriging: the best linear estimate under a fitted variogram, with the variance of its error. The bands below are estimate ± 1.96σ, and they widen wherever the holes thin out.
    
          
      
            
    • The variogram is fitted by weighted least squares over a fixed grid of ranges, exponential or spherical, so the same data always gives the same model.
      • 
            
    • The kriging system is factored once per variable, with LU and partial pivoting, and each query is one solve.
      • 
            
    • The numerics, the variogram fit and the random number generator are all in the crate. It has no dependencies, and every result is reproducible from a seed.
      • 
          
    
        
    
        
    
          
    
            
    
              
    Section A to B, 944 m through six holes
    
              
    
                
                
                
              
    
            
    
            
    
            
    
            
    95 stations along the line from (40, 470) to (980, 560), each one strata query --json on examples/site.csv, with no σ limit so every band shows. Holes within 60 m are projected onto the line, as strata section does. The site is synthetic, so the truth is known: the dashed lines are the seed-7 geology, from the crate's own Truth::new(7). At all 475 points drawn, the truth is inside the band.
    
          
    
        
    
      
    
    
    

    
    
      
    
        
    
          
    4Measure the variance before trusting it
    
          
    Kriging variance assumes the fitted variogram is exactly right, and a variogram fitted to forty holes isn't. So each variable is measured before it's used: each hole is kriged from the others, and the variance is scaled by the mean squared standardised error.
    
          
      
            
    • On smooth geology the fitted models overstate short-range variation, so raw intervals are too wide. On rough geology the sill comes out a little low, so they are slightly narrow.
      • 
            
    • The factor is clamped to [0.25, 4], so a handful of holes can't make the model wildly sure or wildly unsure.
      • 
            
    • strata validate prints every factor and marks any that hit a limit. On the example site, the ground elevation does.
      • 
          
    
        
    
        
    
          
    
            
    
              
    100 sites each, before and after the correction
    
              
    
                
                
                
              
    
            
    
            
    
              
    
              
    
                
    
                  
    The query 120 m from a hole
    
                  
    
                
    
                
    
                  
    What strata validate measured on the example
    
                  
    
                
    
              
    
            
    
            
    Measured correction is the current code. Raw is the same code with the factor fixed at 1, built in a scratch copy only for this comparison. Both columns are strata bench --sites 100, smooth and --rough; the shaded band is the 92 to 98% that CI's --check accepts. The validate table is verbatim; the asterisk marks a factor at its limit.
    
          
    
        
    
      
    
    
    

    
    
      
    
        
    
          
    5Refuse where the interval is too wide
    
          
    This is the line. A contact is answered only if its σ is at most --max-sigma, 1.5 m by default, which is a 95% interval of about ± 2.9 m. Past that, strata says REFUSED and gives the σ instead of a number someone might design to.
    
          
      
            
    • Errors only accumulate downward, so once one contact is refused, every contact below it is refused too. An answer is always a prefix of the column.
      • 
            
    • A unit that fewer than six holes got through isn't modelled at all, and nothing below it is answered.
      • 
            
    • A looser limit never answers less. The property tests check that at seven limits, from 0.1 m to no limit at all.
      • 
          
    
          
    The answer to “where do we need another borehole?” is the hatching. The limit is a flag, so the person who carries the risk sets it, not the tool.
    
        
    
        
    
          
    
            
    
              
    The same section under three limits
    
              
    
                
                
                
                
              
    
            
    
            
    
            
    
            
    
            
    The hatching applies the rule in Column::answer to the σ the CLI returned at each station: hatched below the deepest contact it would still answer. The query is real at each limit, at (600, 480), which is 120.7 m from BH16 and 44 m off the section line.
    
          
    
        
    
      
    
    
    

    
    
      
    ClearedAbove this line strata builds and measures the model. Below it, an elevation leaves the tool only with its interval, and the intervals are held to account against holes it never saw.
    
    
    
    
  
    

  
    
    
    
      
    
        
    
          
    6Answer with an exit code a script can't misread
    
          
    strata query exits 1 if anything was refused, so a script can't mistake a partial answer for a full one. Bad input or usage exits 2. The same answers come out as a table, as JSON, or as SVG sections and maps.
    
          
      
            
    • The JSON gives each surface a status: known with its elevation, σ and 95% bounds, too_uncertain with the σ and the limit, or no_data.
      • 
            
    • strata section draws the layers and bands with the refused column hatched; strata map shades σ across the site and hatches where the chosen surface is refused.
      • 
            
    • Rust ignores SIGPIPE, so piping a section into head would panic. strata restores the default and exits quietly, and a test holds it to that.
      • 
          
    
        
    
        
    
          
    
            
    
              
    Four runs, three exit codes
    
              
            
    
            
    
              
    
              
    
            
    
            
    Real runs on examples/site.csv, output verbatim; the exit code is the shell's. Choose a card to see its run.
    
          
    
        
    
      
    
    
    

    
    
      
    
        
    
          
    7Check the intervals against the truth
    
          
    Real boreholes can't show what lies between them, so the intervals are tested two ways. strata bench builds synthetic sites whose geology is known everywhere and scores the model at 200 random points per site, between the holes. strata cv takes each hole out in turn, refits everything without it, and predicts it.
    
          
      
            
    • The truth isn't from the model's family. Smooth sites are Matérn 3/2, which strata can't fit; rough sites are exponential. Both must pass.
      • 
            
    • In the test suite every surface must pass on its own, not only the pooled figure. CI also runs both 100-site benches on every push and fails if the pooled coverage leaves 92 to 98%.
      • 
            
    • Rough contacts are refused far more often, 95% at the base of gravel against 32%. They really are less predictable from 40 holes, and the model says so.
      • 
          
    
          
    Forty holes give a noisy coverage figure: one miss moves it by 2.5 points. So cv --check only fails outside 85 to 99.5%, which still catches a model that's badly off.
    
        
    
        
    
          
    
            
    
              
    Inside the 95% interval, surface by surface
    
              
    
                
                
                
              
    
            
    
            
    
              
    
              
    
            
    
            
    Every number is from the command shown, run on the release binary: strata cv examples/site.csv, and strata bench --sites 100 with and without --rough, 20,000 scored points per surface. The shaded band is what that command's --check accepts. Refused is the share of points whose σ was over 1.5 m; they are still scored.
    
          
    
        
    
      
    
    
    

    
    
      
    
        
    
          
    46 tests, in a few seconds
    
          
    cargo test runs unit tests, properties over random sites, calibration against synthetic truth, and the real binary as a program. CI runs them on Rust 1.74 and on stable.
    
        
    
        
    
          
    
            
            
      
          
      
          
      
            
      Weighted toward the claim
      
            
        
              
      • Properties, on up to 25 random sites: layers never cross anywhere, every hole is honoured exactly, variance only grows downward, a looser limit never answers less, and far from every hole the model refuses.
        • 
              
      • Calibration, 40 sites × 150 points: every surface inside 90 to 98.5%, pooled 93 to 97% on smooth sites and 92 to 97% on rough ones, msse 0.75 to 1.25. More holes must give lower error and fewer refusals.
        • 
              
      • Numerics: LU against random systems, including the zero diagonal of a kriging matrix; kriging exact at the data, weights summing to 1, variance growing with distance.
        • 
              
      • The parser's 19 ways to reject a file, each with its line number.
        • 
              
      • The CLI: exit codes and output, a closed pipe, and a check that examples/site.csv is exactly what strata synth --seed 7 produces.
        • 
            
      
            
      CI also runs rustfmt and clippy with warnings as errors, both 100-site benches, and a check that examples/ and docs/ match what the current code generates.
      
          
      
        
      
      
      
    
      

    
      
      
      
        
      
          
      
            
      Deliberately missing
      
            
      
              
      Faults, lenses and repeated units
      Layer-cake stratigraphy only: each unit is one layer, in one order. The parser rejects a unit logged twice in a hole instead of guessing.
      
              
      Cokriging
      Thicknesses are kriged independently, and adding their variances assumes uncorrelated errors. The leave-one-out correction and strata cv catch the effect on the intervals; cokriging would model it.
      
              
      Non-stationary variograms
      One variogram per variable for the whole site. A site that is really two settings should be split into two files.
      
              
      Censored data
      A hole that stopped early contributes nothing to the unit it stopped in. It is checked against the model, not used to fit it.
      
              
      A non-Gaussian 95%
      Intervals are estimate ± 1.96σ. The benches show this holds for the fields tested, not for every real site, so run strata cv on yours.
      
            
      
          
      
          
      
            
      Run it
      
            
      git clone https://github.com/ahmed-hashim-pro/strata && cd strata
cargo build --release
export PATH=$PWD/target/release:$PATH

strata validate examples/site.csv
strata query    examples/site.csv 490 440
strata section  examples/site.csv --from 40,470 --to 980,560 -o section.svg
strata cv       examples/site.csv
strata bench

# The full suite, then the checks CI adds.
cargo test
make check
      
            
      Rust 1.74 or newer and nothing else: the crate has no dependencies. Depths are vertical, so inclined holes need converting to true vertical depth first.