### Visualize Geospatial Data with viz

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Visualizes geospatial data using the `viz` recipe function from GeoStats.jl. This example visualizes the geometry of `Ω` and colors it based on the associated `z` data.

```julia
viz(Ω.geometry, color = Ω.z)
```

--------------------------------

### GeoStats.jl Interpolation Setup

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/faq/interpolation.md

Initializes the GeoStats and CairoMakie libraries for interpolation tasks. This is a common setup step for performing spatial interpolation with GeoStats.jl.

```julia
using GeoStats # hide
import CairoMakie as Mke # hide
```

--------------------------------

### View Geospatial Data with viewer

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Launches an interactive scientific viewer for geospatial data using the `viewer` function. This viewer displays all viewable variables with a colorbar and other interactive controls.

```julia
viewer(Ω)
```

--------------------------------

### Geostatistical Interpolation Example with GeoStats.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/index.md

This example demonstrates a basic geostatistical interpolation using Kriging on a Cartesian grid. It involves loading the framework, defining data, choosing an interpolation model, performing the interpolation, and visualizing the results.

```julia
# load framework
using GeoStats

# load visualization backend
import CairoMakie as Mke

# attribute table
table = (; z=[1.,0.,1.])

# coordinates for each row
coord = [(25.,25.), (50.,75.), (75.,50.)]

# georeference data
geotable = georef(table, coord)

# interpolation domain
grid = CartesianGrid(100, 100)

# choose an interpolation model
model = Kriging(GaussianVariogram(range=35.))

# perform interpolation over grid
interp = geotable |> Interpolate(grid, model=model)

# visualize the solution
interp |> viewer
```

--------------------------------

### Install GeoStats.jl Package Manager

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/index.md

This snippet shows how to install the latest stable release of the GeoStats.jl package using Julia's built-in package manager.

```julia
using Pkg
Pkg.add("GeoStats")
```

--------------------------------

### Calculate Area and Perimeter from Meshes.jl Geometries

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Demonstrates how to calculate the total area and perimeter of geometries loaded using Meshes.jl. These functions are part of the Meshes.jl package and operate on the geometry objects.

```julia
sum(area, zone.geometry)
```

```julia
sum(perimeter, zone.geometry)
```

--------------------------------

### Geospatial Queries with @transform Macro in Geostats.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Shows an example of using the @transform macro from Geostats.jl for geospatial data manipulation. This macro simplifies operations involving the 'geometry' column, such as calculating area.

```julia
@transform(Ω, :w = 2 * :z * area(:geometry))

```

--------------------------------

### Create and Visualize a Geotable for Clustering

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/clustering.md

This example demonstrates the creation of a sample geotable with synthetic data and visualizes it. This geotable is used to illustrate the geostatistical clustering methods provided by the package.

```julia
gtb = georef((z=[10sin(i/10) + j for i in 1:4:100, j in 1:4:100],))

gt b |> viewer
```

--------------------------------

### GeoTable Creation

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/faq/coords.md

Creates a sample GeoTable with random data for demonstration purposes. This is a common starting point for applying transformations.

```julia
data = georef((; z=rand(1000, 1000)))
```

--------------------------------

### Accessing Table Data in Geostats.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Demonstrates accessing data from a geospatial dataset that conforms to the Tables.jl interface. It shows how to access specific rows and columns, including a dynamically generated 'geometry' column.

```julia
Ω[1,:]

```

```julia
Ω[1,"geometry"]

```

```julia
Ω.z

```

```julia
Ω.geometry

```

```julia
values(Ω)

```

```julia
domain(Ω)

```

--------------------------------

### Get GeoStats.jl Version Information (Julia)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/contributing.md

These commands retrieve the version information for Julia and the installed packages, including GeoStats.jl, to aid in reporting issues.

```julia
julia> versioninfo()
julia> using Pkg; Pkg.status()
```

--------------------------------

### Example Partitioning with GeoStats.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/partitioning.md

Demonstrates how to use the `partition` function with `BisectFractionPartition` to split a geotable into two parts based on a specified direction and fraction. It also shows how to access the partitioned results and use the `geosplit` utility function as a shortcut.

```julia
using GeoStats # hide
import CairoMakie as Mke # hide

# random image with 100x100 pixels
gtb = georef((; z=rand(100, 100)))

# 70%/30% partition perpendicular to (1.0, 1.0) direction
Π = partition(gtb, BisectFractionPartition((1.0, 1.0), fraction=0.7))

Π[1] |> viewer

Π[2] |> viewer

geosplit(gtb, 0.7, (1.0, 1.0))
```

--------------------------------

### Applying Geospatial Transforms with Pipelines in Geostats.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Illustrates how to create and apply geospatial transform pipelines in Geostats.jl. It shows combining transforms, feeding data to the pipeline, and visualizing the data distribution before and after the transform. It also demonstrates checking the bounding box of geometries before and after.

```julia
# pipeline with table transforms
pipe = Quantile() → StdCoords()

# feed geospatial data to pipeline
Ω̂ = pipe(Ω)

# plot distribution before and after pipeline
fig = Mke.Figure(size = (800, 400))
Mke.hist(fig[1,1], Ω.z, color = :gray)
Mke.hist(fig[2,1], Ω̂.z, color = :gray)
fig

```

```julia
boundingbox(Ω.geometry)

```

```julia
boundingbox(Ω̂.geometry)

```

--------------------------------

### Initialize Geostats and Plotting - Julia

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/interpolation.md

This code snippet initializes the GeoStats package and imports CairoMakie for plotting. It's a common setup for geostatistical visualization tasks in Julia.

```julia
using GeoStats
import CairoMakie as Mke
```

--------------------------------

### CoKriging Multivariate Interpolation (Heterotopic)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/interpolation.md

This example shows how to perform CoKriging in a heterotopic setting, where missing values indicate different spatial supports. It uses a similar setup to the homotopic example but includes missing values in the data table, which are automatically detected by the interpolation process.

```julia
# 3 variables with missing values
table = (; a=[1.0, 0.0, 0.0], b=[0.0, 1.0, missing], c=[missing, 0.0, 1.0])

gtb = georef(table, coord)

itp = gtb |> Interpolate(grid, model=Kriging(γ))
itp |> Select("a") |> viewer
itp |> Select("b") |> viewer
itp |> Select("c") |> viewer
```

--------------------------------

### Define and Train DecisionTreeClassifier Model

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Defines features and the label for a geostatistical learning model and initializes a DecisionTreeClassifier. This model is then used with the `Learn` transform to fit the training data (Ωs). Any model from StatsLearnModels.jl or ScikitLearn.jl is compatible.

```julia
feats = ["band1", "band2", "band3", "band4"]
label = "crop"

model = DecisionTreeClassifier()

learn = Learn(Ωs, model, feats => label)
```

--------------------------------

### Visualize Model Predictions Side-by-Side with True Labels

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Visualizes the model's predictions (`Ω̂t`) and the true labels (`Ωt.crop`) side-by-side for comparison. This allows for a direct assessment of the model's performance in predicting the crop type across different locations. Requires Makie.jl.

```julia
fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], Ω̂t.geometry, color = Ω̂t.crop)
viz(fig[1,2], Ωt.geometry, color = Ωt.crop)
fig
```

--------------------------------

### Load and Inspect Geospatial Data with CSV.jl and DataFrames.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Loads geospatial data from a CSV file into a DataFrame and displays the first 5 rows. This step is crucial for understanding the structure and content of the input data before georeferencing and modeling. It requires the CSV.jl and DataFrames.jl packages.

```julia
using CSV
using DataFrames

table = CSV.File("data/agriculture.csv") |> DataFrame

first(table, 5)
```

--------------------------------

### Georeference and Visualize Geospatial Data with Geostats.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Georeferences a DataFrame by specifying 'x' and 'y' columns as coordinates and then visualizes two variables ('band4' and 'crop') on separate plots. This prepares the data for geospatial analysis and visual inspection. Requires Makie.jl (aliased as Mke) and Geostats.jl.

```julia
Ω = georef(table, ("x", "y"))

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], Ω.geometry, color = Ω.band4)
viz(fig[1,2], Ω.geometry, color = Ω.crop)
fig
```

--------------------------------

### Simulating a Gaussian Process Field

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/simulation.md

This example demonstrates how to simulate a Gaussian Process, a type of field process in GeoStats.jl. It defines a grid, a GaussianProcess with a specified variogram, and then generates 100 realizations of the process over the grid.

```julia
# domain of simulation
grid = CartesianGrid(100, 100)

# field process
proc = GaussianProcess(GaussianVariogram(range=30.0))

# perform simulation
real = rand(proc, grid, 100)
```

--------------------------------

### Auxiliary Variable Simulation - GeoStats.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/fieldprocs.md

Demonstrates incorporating auxiliary variables to guide pattern selection in spatial simulations. This example uses a blur filter as a forward model to generate auxiliary data. Requires ImageFiltering.jl, GeoStats.jl, and GeoStatsImages.jl.

```julia
using ImageFiltering
using GeoStats
using GeoStatsImages

# image assumed as ground truth (unknown)
truth = geostatsimage("WalkerLakeTruth")

# training image with similar patterns
img = geostatsimage("WalkerLake")

# forward model (blur filter)
function forward(data)
  img = asarray(data, :Z)
  krn = KernelFactors.IIRGaussian([10,10])
  fwd = imfilter(img, krn)
  georef((; fwd=vec(fwd)), domain(data))
end
```

--------------------------------

### Create Georeferenced Data with georef

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Creates georeferenced geospatial data from a Julia array. The `georef` function assigns default coordinates based on the array dimensions, creating an object `Ω` with associated data `z`.

```julia
z = [10sin(i/10) + 2j for i in 1:50, j in 1:50]

Ω = georef((; z=z))
```

--------------------------------

### Calculating and Visualizing Ensemble Statistics

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/simulation.md

This example shows how to calculate and visualize the mean and variance of an ensemble of simulated field processes. It computes the mean and variance of the 'real' simulation results and then visualizes them using CairoMakie.

```julia
m, v = mean(real), var(real)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], m.geometry, color = m.field)
viz(fig[1,2], v.geometry, color = v.field)
fig
```

--------------------------------

### Generate Predictions with Trained Model

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Generates predictions on the test set (Ωt) using the previously trained `Learn` transform. The output `Ω̂t` is a geospatial data object containing the predictions. This step applies the learned model to new, unseen data.

```julia
Ω̂t = learn(Ωt)
```

--------------------------------

### Creating HScatter Plots in GeoStats.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/visualization.md

Provides an example of creating multiple hscatter plots using GeoStats.jl and CairoMakie.jl. Hscatter plots visualize the relationship between two variables at different lags.

```julia
using GeoStats # hide
import CairoMakie as Mke # hide

𝒟 = georef((Z=[10sin(i/10) + j for i in 1:100, j in 1:200],))

𝒮 = sample(𝒟, UniformSampling(500))

fig = Mke.Figure(size = (800, 400))
hscatter(fig[1,1], 𝒮, :Z, :Z, lag=0)
hscatter(fig[1,2], 𝒮, :Z, :Z, lag=20)
hscatter(fig[2,1], 𝒮, :Z, :Z, lag=40)
hscatter(fig[2,2], 𝒮, :Z, :Z, lag=60)
fig
```

--------------------------------

### Thinning a Binomial Point Pattern (Julia)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/simulation.md

Demonstrates thinning a pre-existing Binomial point pattern with a specified probability using `thin` and `RandomThinning`. The example shows the original pattern and the thinned version, both visualized within a box geometry.

```julia
# geometry of interest
box = Box((0, 0), (100, 100))

# Binomial process
proc = BinomialProcess(2000)

# sample point pattern
pset₁ = rand(proc, box)

# thin point pattern with probability 0.5
pset₂ = thin(pset₁, RandomThinning(0.5))

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], box)
viz!(fig[1,1], pset₁, color = :black)
viz(fig[1,2], box)
viz!(fig[1,2], pset₂, color = :black)
fig
```

--------------------------------

### Visualize Geospatial Data

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/functions.md

Generates a 2D Cartesian grid of points with associated random Z values and visualizes them using `viz`. This sets up data for interpolation examples.

```julia
using Random # hide
Random.seed!(2000) # hide

points = rand(Point, 50, crs=Cartesian2D) |> Scale(100)
geotable = georef((Z=rand(50),), points)

viz(geotable.geometry, color = geotable.Z)
```

--------------------------------

### Work with Remote Datasets and Visualize Layers using GeoArtifacts

Source: https://context7.com/juliaearth/geostats.jl/llms.txt

This snippet illustrates how to access and visualize remote geospatial datasets using the GeoArtifacts.jl package. It includes examples for loading natural earth data and visualizing multiple layers by chaining transformations and using the viewer function. Dependencies include GeoStats.jl, GeoIO.jl, and GeoArtifacts.jl.

```julia
using GeoStats
using GeoIO
using GeoArtifacts

# Work with remote datasets
earth = NaturalEarth.naturalearth1("water")
borders = NaturalEarth.borders()
airports = NaturalEarth.airports()

# Visualize multiple layers
earth |> Upscale(10, 5) |> viewer
viz!(borders.geometry, color="cyan")
viz!(airports.geometry, color="black", pointsize=4)
```

--------------------------------

### Split Geospatial Data with geosplit Function

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/quickstart.md

Splits the georeferenced data (Ω) into training (Ωs) and testing (Ωt) sets using the `geosplit` function. A separating plane is defined by its normal direction (1.0, -1.0), and 20% of the data is reserved for the training set. The function also visualizes the domains of both sets. Requires Makie.jl and Geostats.jl.

```julia
Ωs, Ωt = geosplit(Ω, 0.2, (1.0, -1.0))

viz(Ωs.geometry)
viz!(Ωt.geometry, color = :gray90)
Mke.current_figure()
```

--------------------------------

### Load Geospatial Data from File (Shapefile)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/data.md

Loads geospatial data from various file formats using the GeoIO.jl package. This example specifically loads a shapefile. The loaded data can then be used with Geostats.jl functions.

```julia
using GeoIO

zone = GeoIO.load("data/zone.shp")
```

--------------------------------

### Download Geospatial Data Artifacts

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/data.md

Downloads geospatial data from online repositories using the GeoArtifacts.jl package. This example downloads natural earth data like water bodies, borders, airports, and ports, and visualizes them.

```julia
using GeoArtifacts
import CairoMakie as Mke

earth    = NaturalEarth.naturalearth1("water")
borders  = NaturalEarth.borders()
airports = NaturalEarth.airports()
ports    = NaturalEarth.ports()

earth |> Upscale(10, 5) |> viewer

viz!(borders.geometry, color = "cyan")
viz!(airports.geometry, color = "black", pointsize=4, pointmarker='✈')
viz!(ports.geometry, color = "blue", pointsize=4, pointmarker='⚓')

Mke.current_figure()
```

--------------------------------

### Visualizing Cumulative Distribution Functions

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/simulation.md

This example illustrates how to compute and visualize the cumulative distribution functions (CDFs) at specific percentiles for an ensemble of simulated field processes. It calculates the CDFs at 0.25 and 0.75 and visualizes them.

```julia
p25 = cdf(real, 0.25)
p75 = cdf(real, 0.75)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], p25.geometry, color = p25.field)
viz(fig[1,2], p75.geometry, color = p75.field)
fig
```

--------------------------------

### Multivariate Gaussian Process Simulation

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/simulation.md

This example demonstrates multivariate geostatistical simulation (cosimulation) with the GaussianProcess. It defines a multivariate covariance function to model the correlation between variables and then simulates two variables jointly, visualizing them side-by-side.

```julia
# multivariate covariance function
func = [1.0 0.8; 0.8 1.0] * SphericalCovariance(range=30.0)

# multivariate Gaussian process
proc = GaussianProcess(func)

# simulate 2 variables jointly
real = rand(proc, grid)

# visualize variables side by side
fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real.geometry, color = real.field1)
viz(fig[1,2], real.geometry, color = real.field2)
fig
```

--------------------------------

### Combine Point Processes with Union (Julia)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/simulation.md

Demonstrates the superposition of two Binomial point processes using the union operator (∪). This operation combines points from both processes. The input is two BinomialProcess objects, and the output is a new BinomialProcess representing their union. The example visualizes the resulting point patterns on a box geometry.

```julia
# geometry of interest
box = Box((0, 0), (100, 100))

# superposition of two Binomial processes
proc₁ = BinomialProcess(500)
proc₂ = BinomialProcess(500)
proc  = proc₁ ∪ proc₂ # 1000 points

pset = rand(proc, box, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], box)
viz!(fig[1,1], pset[1], color = :black)
viz(fig[1,2], box)
viz!(fig[1,2], pset[2], color = :black)
fig
```

--------------------------------

### Thinning a Poisson Process (Julia)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/simulation.md

Illustrates how to reduce the intensity of a Poisson point process by half using the `thin` function with `RandomThinning`. The example compares a thinned process with the original one and visualizes their sampled point patterns.

```julia
# geometry of interest
box = Box((0, 0), (100, 100))

# reduce intensity of Poisson process by half
proc₁ = PoissonProcess(0.5)
proc₂ = thin(proc₁, RandomThinning(0.5))

# sample point patterns
pset₁ = rand(proc₁, box)
pset₂ = rand(proc₂, box)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], box)
viz!(fig[1,1], pset₁, color = :black)
viz(fig[1,2], box)
viz!(fig[1,2], pset₂, color = :black)
fig
```

--------------------------------

### Universal Kriging with Trend in Julia

Source: https://context7.com/juliaearth/geostats.jl/llms.txt

Demonstrates Universal Kriging, which models spatial trends using polynomial drift functions. The examples show how to set up Universal Kriging models with linear (degree=1) and quadratic (degree=2) trends using the GeoStats library.

```julia
using GeoStats

# Sample data with spatial trend
x_coords = 0:0.5:10
y_coords = 0:0.5:10
trend = [x + 0.5*y + 0.1*x*y for x in x_coords, y in y_coords]
noise = 0.1 * randn(length(x_coords), length(y_coords))
z_values = trend + noise

data = georef((Z=vec(z_values),))

# Define Universal Kriging with polynomial degree
uk_degree1 = Kriging(
    GaussianVariogram(range=2.0),
    degree=1  # Linear trend
)

uk_degree2 = Kriging(
    GaussianVariogram(range=2.0),
    degree=2  # Quadratic trend
)

```

--------------------------------

### Inhibition Process Sampling in Julia

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/pointprocs.md

Provides an example of sampling from an InhibitionProcess in Julia using GeoStats.jl. This process model ensures a minimum distance between points. The code defines the process and generates sample point patterns within a specified box, visualized using CairoMakie.

```julia
using GeoStats # hide
import CairoMakie as Mke # hide

# geometry of interest
box = Box((0, 0), (100, 100))

# inhibition process
proc = InhibitionProcess(2.0)

# sample point pattern
pset = rand(proc, box, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], box)
viz!(fig[1,1], pset[1], color = :black)
viz(fig[1,2], box)
viz!(fig[1,2], pset[2], color = :black)
fig
```

--------------------------------

### Apply GHC Clustering to Geotable

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/clustering.md

This example demonstrates the application of the GHC (Geostatistical Hierarchical Clustering) algorithm to a geotable. It clusters the data into 5 regions with a smoothing factor of 1.0, considering both values and domain, and visualizes the result.

```julia
ctb = gtb |> GHC(5, 1.0)

ctb |> viewer
```

--------------------------------

### CoKriging Multivariate Interpolation (Homotopic)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/interpolation.md

This example demonstrates CoKriging for multivariate interpolation with three variables. It defines a table of data, coordinates, and a multivariate variogram model. The code then performs interpolation on a grid and selects individual variables for visualization.

```julia
# 3 variables
table = (; a=[1.0, 0.0, 0.0], b=[0.0, 1.0, 0.0], c=[0.0, 0.0, 1.0])
coord = [(25.0, 25.0), (50.0, 75.0), (75.0, 50.0)]

# 3x3 variogram
γ = [1.0 0.3 0.1; 0.3 1.0 0.2; 0.1 0.2 1.0] * SphericalVariogram(range=35.)

gtb = georef(table, coord)

itp = gtb |> Interpolate(grid, model=Kriging(γ))
itp |> Select("a") |> viewer
itp |> Select("b") |> viewer
itp |> Select("c") |> viewer
```

--------------------------------

### Apply SLIC Clustering to Geotable

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/clustering.md

This example applies the SLIC (Simple Linear Iterative Clustering) algorithm to a geotable. It divides the data into 5 clusters with a compactness factor of 0.01, considering spatial proximity, and visualizes the clustered output.

```julia
ctb = gtb |> SLIC(5, 0.01)

ctb |> viewer
```

--------------------------------

### Universal Kriging Mathematical Formulation

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/interpolation.md

This section details the mathematical setup for Universal Kriging, including assumptions about the random field mean, polynomial matrices, variogram matrices, and the resulting linear system. It also defines how the mean and variance at a specific location are calculated.

```julia
In Universal Kriging, the mean of the random field is assumed to be a linear combination of known smooth functions.
For example, it is common to assume

μ(ρ) = ∑_{k=1}^{N_d} β_k f_k(ρ)

with ``N_d`` monomials ``f_k`` of degree up to ``d``. For example, in 2D there are ``6`` monomials of degree up to ``2``:

μ(x_1,x_2) =  β_1 1 + β_2 x_1 + γ x_2 + β_4 x_1 x_2 + β_5 x_1^2 + β_6 x_2^2

The choice of the degree ``d`` determines the size of the polynomial matrix

ℂ =
∣∣∣∣
 f_1(ρ_1) & f_2(ρ_1) & ⋯ & f_{N_d}(ρ_1) \\
f_1(ρ_2) & f_2(ρ_2) & ⋯ & f_{N_d}(ρ_2) \\
⋰ & ⋰ & ddots & ⋰ \\
f_1(ρ_n) & f_2(ρ_n) & ⋯ & f_{N_d}(ρ_n) \\
∣∣∣∣

and polynomial vector ``→ = ∣∣∣
 f_1(ρ_0) & f_2(ρ_0) & ⋯ & f_{N_d}(ρ_0) \\∣∣∣
⊥.

The variogram determines the variogram matrix:

℉ =
∣∣∣∣
 γ(ρ_1,ρ_1) & γ(ρ_1,ρ_2) & ⋯ & γ(ρ_1,ρ_n) \\
 γ(ρ_2,ρ_1) & γ(ρ_2,ρ_2) & ⋯ & γ(ρ_2,ρ_n) \\
 ⋰ & ⋰ & ddots & ⋰ \\
 γ(ρ_n,ρ_1) & γ(ρ_n,ρ_2) & ⋯ & γ(ρ_n,ρ_n) \\
∣∣∣∣

and the variogram vector
``→ = ∣∣∣
 γ(ρ_1,ρ_0) & γ(ρ_2,ρ_0) & ⋯ & γ(ρ_n,ρ_0) \\
∣∣∣
⊥.

The resulting linear system is:

∣∣∣
 ℉ & ℂ \\
 ℂ⁺ & ⊛ \\
∣∣∣
 ∣∣∣
 → \\
 → \\
∣∣∣

with ``∣∣∣ ∂`` the Lagrange multipliers associated with the universal constraints. The mean and
variance at location ``ρ_0`` are given by:

μ(ρ_0) = ⊥⁺ →

σ^2(ρ_0) = ∣∣∣
→ \\
 → \\
∣∣∣⁺ ∣∣∣
 → \\
 ∂ \\
∣∣∣

```

--------------------------------

### Apply Potrace Transform (Julia)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/transforms.md

This example demonstrates using the Potrace transform to trace polygons from a binary mask derived from continuous feature data. It generates continuous data, creates a binary mask, georeferences it, applies Potrace, and visualizes the original data and the traced polygons.

```julia
using GeoStats
using CairoMakie as Mke

# continuous feature
Z = [sin(i/10) + sin(j/10) for i in 1:100, j in 1:100]

# binary mask
M = Z .> 0

# georeference data
Ω = georef((Z=Z, M=M))

# trace polygons using mask
𝒯 = Ω |> Potrace(:M)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], Ω.geometry, color = Ω.Z)
viz(fig[1,2], 𝒯.geometry, color = 𝒯.Z)
fig

𝒯.geometry
```

--------------------------------

### Apply and Revert Quantile Transform (Julia)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/transforms.md

This example shows how to apply a Quantile transform to geospatial data and then revert it. It first georeferences a DataFrame onto a Cartesian grid, applies the Quantile transform to normalize the features, and then demonstrates reverting the transform to recover the original data.

```julia
using DataFrames

# table of features and domain
tab = DataFrame(a=rand(1000), b=randn(1000), c=rand(1000))
dom = CartesianGrid(100, 100)

# georeference table onto domain
Ω = georef(tab, dom)

# describe features
describe(values(Ω))

pipe = Quantile()

Ω̄, cache = apply(pipe, Ω)

describe(values(Ω̄))

Ωₒ = revert(pipe, Ω̄, cache)

describe(values(Ωₒ))
```

--------------------------------

### Plot Surfaces of Association

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/functions.md

Documents the `surfplot` and `surfplot!` functions for plotting surfaces of association given a normal direction. An example is provided for an anisotropic variogram.

```julia
surfplot
surfplot!

# example
γ = GaussianVariogram(ranges=(3, 2, 1))
surfplot(γ)
```

--------------------------------

### Get Raster Grid Size

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/faq/rasters.md

Retrieves and displays the dimensions of the underlying grid for a given raster. This is useful for understanding the spatial resolution and extent of the data.

```julia
size(raster.geometry)
```

--------------------------------

### Plot Geostatistical Functions

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/functions.md

Provides documentation for the `funplot` and `funplot!` functions, which are used to plot geostatistical functions. An example is given for plotting an anisotropic Gaussian variogram.

```julia
funplot
funplot!

# example
γ = GaussianVariogram(ranges=(3, 2, 1))
funplot(γ)
```

--------------------------------

### GeoTable Translation Transform

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/faq/coords.md

Applies a `Translate` transformation to a GeoTable, shifting its geometries by specified distances in meters and feet. This is an example of rigid motion.

```julia
data |> Translate(10u"m", 8u"ft")
```

--------------------------------

### Visualizing Gaussian Process Realizations

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/simulation.md

This code snippet visualizes the first two realizations of a simulated Gaussian Process. It uses CairoMakie to create a figure and display the simulated field values on their respective geometries.

```julia
fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].field)
viz(fig[1,2], real[2].geometry, color = real[2].field)
fig
```

--------------------------------

### Compare Kriging and Inverse Distance Weighting Solvers in GeoStats.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/paper/paper.md

This snippet illustrates how to compare different geostatistical solvers, specifically Kriging and Inverse Distance Weighting, on a predefined estimation problem. It utilizes the `compare` function with `VisualComparison` to visually assess the performance of each solver. Dependencies include GeoStats.jl and Plots.jl. The input is a problem definition and a list of solvers; the output is a visual comparison of their solutions.

```julia
using GeoStats
using Plots

# define solvers to be compared
solver1 = Kriging(
    :precipitation => @NT(variogram=GaussianVariogram(range=35.))
)
solver2 = InvDistWeight()

# compare solvers with a comparison method (e.g. visual comparison)
compare([solver1, solver2], problem, VisualComparison())
```

--------------------------------

### GeoTable Projection Transform

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/faq/coords.md

Reprojects a GeoTable into a new coordinate reference system (CRS) using the `Proj` function. This example uses `Polar` as the target CRS.

```julia
data |> Proj(Polar)
```

--------------------------------

### Binomial Process Sampling in Julia

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/pointprocs.md

Demonstrates how to create and sample from a BinomialProcess in Julia using the GeoStats.jl package. This involves defining a bounding box and then generating point patterns within that box. It utilizes CairoMakie for visualization.

```julia
using GeoStats # hide
import CairoMakie as Mke # hide

# geometry of interest
box = Box((0, 0), (100, 100))

# Binomial process
proc = BinomialProcess(1000)

# sample point patterns
pset = rand(proc, box, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], box)
viz!(fig[1,1], pset[1], color = :black)
viz(fig[1,2], box)
viz!(fig[1,2], pset[2], color = :black)
fig
```

--------------------------------

### Stratigraphic Simulation with StratiGraphics.jl (Julia)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/fieldprocs.md

Simulates stratigraphic environments using StratiGraphics.jl. It sets up a grid, defines a smoothing process and environmental parameters, and then generates stratigraphic simulations, visualizing their geometry and field values.

```julia
using StratiGraphics

# domain of interest
grid = CartesianGrid(50, 50, 20)

# stratigraphic environment
p = SmoothingProcess()
T = [0.5 0.5; 0.5 0.5]
Δ = ExponentialDuration(1.0)
ℰ = Environment([p, p], T, Δ)

# strata simulation
real = rand(StrataProcess(ℰ), grid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].field)
viz(fig[1,2], real[2].geometry, color = real[2].field)
fig
```

--------------------------------

### Fit Theoretical Functions to Empirical Data

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/functions.md

Documents the `fit` function from GeoStatsFunctions, which is used to fit theoretical geostatistical functions to empirical estimates. An example of creating sinusoidal data for an empirical variogram is provided.

```julia
GeoStatsFunctions.fit

# example
# sinusoidal data
𝒟 = georef((Z=[sin(i/2) + sin(j/2) for i in 1:50, j in 1:50],))
```

--------------------------------

### Create and Visualize PointSet

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/domains.md

Creates a PointSet domain with 100 random points and visualizes it. A PointSet is used to represent discrete locations like rain gauge stations.

```julia
pset = PointSet(rand(Point, 100))

viz(pset)
```

--------------------------------

### Quilting Process Simulation - ImageQuilting.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/fieldprocs.md

Simulates unconditional and conditional realizations using the QuiltingProcess from ImageQuilting.jl. This process generates spatial fields based on patterns from a training image. It supports simulation on grids, views of grids, and conditioning on data. Requires ImageQuilting.jl, GeoStats.jl, and GeoStatsImages.jl.

```julia
using ImageQuilting
using GeoStatsImages
using GeoStats
import CairoMakie as Mke

# domain of interest
grid = CartesianGrid(200, 200)

# quilting process
img  = geostatsimage("Strebelle")
proc = QuiltingProcess(img, (62, 62))

# unconditional simulation
real = rand(proc, grid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].facies)
viz(fig[1,2], real[2].geometry, color = real[2].facies)
fig
```

```julia
using ImageQuilting
using GeoStatsImages
using GeoStats
import CairoMakie as Mke

# domain of interest
grid = CartesianGrid(200, 200)

# quilting process
img  = geostatsimage("StoneWall")
proc = QuiltingProcess(img, (13, 13))

# unconditional simulation
real = rand(proc, grid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].Z)
viz(fig[1,2], real[2].geometry, color = real[2].Z)
fig
```

```julia
using ImageQuilting
using GeoStatsImages
using GeoStats
import CairoMakie as Mke

# domain of interest
grid = domain(img) # Assuming img is already defined from previous example

# pre-existing observations
img  = geostatsimage("Strebelle")
data = img |> Sample(20, replace=false)

# quilting process
proc = QuiltingProcess(img, (30, 30))

# conditional simulation
real = rand(proc, grid, 2, data=data)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].facies)
viz(fig[1,2], real[2].geometry, color = real[2].facies)
fig
```

```julia
using ImageQuilting
using GeoStatsImages
using GeoStats
import CairoMakie as Mke

# domain of training image
img = geostatsimage("Strebelle")
grid = domain(img)

# view pixels inside ball
ball = Ball((50.0, 50.0), 25.0)
vgrid = view(grid, ball)

# quilting process
proc = QuiltingProcess(img, (62, 62))

# unconditional simulation
real = rand(proc, vgrid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].facies)
viz(fig[1,2], real[2].geometry, color = real[2].facies)
fig
```

--------------------------------

### Gaussian Process Simulation - GeoStats.jl

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/fieldprocs.md

Simulates unconditional realizations of a Gaussian process on a Cartesian grid using a Gaussian variogram. Requires GeoStats.jl and CairoMakie.jl for visualization. Outputs realizations of the field.

```julia
using GeoStats
import CairoMakie as Mke

# domain of interest
grid = CartesianGrid(100, 100)

# Gaussian process
proc = GaussianProcess(GaussianVariogram(range=30.0))

# unconditional simulation
real = rand(proc, grid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].field)
viz(fig[1,2], real[2].geometry, color = real[2].field)
fig
```

--------------------------------

### Create and Visualize GeometrySet

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/domains.md

Creates a GeometrySet containing a Triangle and a Quadrangle, then visualizes it with segments shown. GeometrySet is suitable for representing collections of GIS geometries.

```julia
tria = Triangle((0.0, 0.0), (1.0, 1.0), (0.0, 1.0))
quad = Quadrangle((1.0, 1.0), (2.0, 1.0), (2.0, 2.0), (1.0, 2.0))
gset = GeometrySet([tria, quad])

viz(gset, showsegments = true)
```

--------------------------------

### Apply Forward Model and Simulate Data (Julia)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/fieldprocs.md

Applies a forward model to images and simulates realizations using QuiltingProcess. This involves generating synthetic data and visualizing the geometric properties and Z values of the simulated realizations.

```julia
# apply forward model to both images
data   = forward(truth)
dataTI = forward(img)

proc = QuiltingProcess(img, (27, 27), soft=(data, dataTI))

real = rand(proc, domain(truth), 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].Z)
viz(fig[1,2], real[2].geometry, color = real[2].Z)
fig
```

--------------------------------

### Build Data Transformation Pipelines with Chain.jl

Source: https://context7.com/juliaearth/geostats.jl/llms.txt

This code demonstrates building complex data transformation pipelines using GeoStats.jl and Chain.jl. It includes steps for filtering outliers, normalizing variables, performing interpolation, and selecting specific variables. It shows both the @chain macro and the arrow operator (→) for pipeline construction. Dependencies include GeoStats.jl and Chain.jl.

```julia
using GeoStats
using Chain

# Load and prepare data
raw_data = georef((
    temperature=rand(100) .* 100,
    pressure=rand(100) .* 1000,
    elevation=rand(100) .* 500
))

# Build transformation pipeline
processed = @chain raw_data begin
    # Remove outliers
    Filter(row -> 10 < row.temperature < 90)

    # Normalize variables
    Normalize()

    # Perform interpolation
    Interpolate(CartesianGrid(50, 50), model=Kriging(GaussianVariogram(range=20.0)))

    # Select specific variables
    Select("temperature", "pressure")

    # Visualize result
    viewer
end

# Alternative pipeline syntax
pipeline = Filter(row -> row.temperature > 0) →
           Normalize() →
           Interpolate(CartesianGrid(50, 50), model=IDW())

result = raw_data |> pipeline
```

--------------------------------

### Kriging Interpolation with GeoStats.jl

Source: https://context7.com/juliaearth/geostats.jl/llms.txt

This snippet illustrates how to perform Kriging interpolation using GeoStats.jl. It involves defining sample data, specifying a target interpolation domain, creating a Kriging model with a variogram, and executing the interpolation. It also shows basic visualization and probabilistic prediction.

```julia
using GeoStats
import CairoMakie as Mke

# Sample data with measurements at specific locations
sample_table = (z=[1.0, 0.0, 1.0])
sample_coords = [(25.0, 25.0), (50.0, 75.0), (75.0, 50.0)]
samples = georef(sample_table, sample_coords)

# Define target interpolation domain
grid = CartesianGrid(100, 100)

# Create Kriging model with Gaussian variogram
model = Kriging(GaussianVariogram(range=35.0))

# Perform interpolation
interpolated = samples |> Interpolate(grid, model=model)

# Visualize results
interpolated |> viewer

# Advanced: probabilistic prediction with variance
using GeoStatsModels
fitted_model = GeoStatsModels.fit(model, samples)
pred_mean, pred_var = GeoStatsModels.predictprob(fitted_model, grid)
```

--------------------------------

### Apply UniqueCoords Transform (Julia)

Source: https://github.com/juliaearth/geostats.jl/blob/master/docs/src/transforms.md

This example demonstrates the use of the UniqueCoords transform to remove duplicate points from a georeferenced dataset. It creates a point set with repeated points, applies the UniqueCoords transform, and shows the resulting dataset with unique coordinates.

```julia
using GeoStats

# point set with repeated points
p = rand(Point, 50)
Ω = georef((Z=rand(100),), [p; p])

# discard repeated points
𝒰 = Ω |> UniqueCoords()
```

--------------------------------

### Neighborhood Search and Interpolation in Julia

Source: https://context7.com/juliaearth/geostats.jl/llms.txt

Demonstrates various neighborhood-based interpolation methods including Nearest Neighbor, Inverse Distance Weighting (IDW), Locally Weighted Regression (LWR), and Polynomial interpolation. It also shows how to perform kriging using a specified neighborhood search. Requires GeoStats and Meshes libraries.

```julia
using GeoStats

# Sample data
table = (value=rand(100),)
points = rand(Point, 100) |> Scale(100)
samples = georef(table, points)
grid = CartesianGrid(50, 50)

# Nearest neighbor interpolation
nn_interp = samples |> Interpolate(grid, model=NN())

# Inverse distance weighting
idw_interp = samples |> Interpolate(grid, model=IDW(power=2.0))

# Locally weighted regression
lwr_interp = samples |> Interpolate(grid, model=LWR(degree=2))

# Polynomial interpolation
poly_interp = samples |> Interpolate(grid, model=Polynomial(degree=2))

# Neighborhood search with kriging
using Meshes
neighborhood = BallNeighborhood(grid, 20.0u"m")  # Search within 20m radius
kriging_model = Kriging(GaussianVariogram(range=15.0))
local_interp = samples |> InterpolateNeighbors(
    grid,
    neighborhood,
    model=kriging_model
)

# Visualize results
local_interp |> viewer
```