GeoStats.jl

GeoStats.jl is an extensible framework for geospatial data science and geostatistical modeling written in Julia. It provides a suite of modules designed for:

• Advanced geometric processing
• State-of-the-art geostatistical algorithms
• Sophisticated visualization of geospatial data

For detailed usage, API references, and guides, refer to the official online documentation.

The Julia geospatial ecosystem is organized into three primary initiatives, each addressing different modeling challenges:

1. JuliaEarth: Focused on geospatial data science and geostatistical modeling. It is designed to be 100% Julia, meaning it avoids reliance on external C/C++ libraries like GDAL or Proj4. It is optimized for both 2D and 3D geometries.
2. JuliaClimate: Focused on climate modeling and providing access to climate data. Key packages include INMET.jl and CDSAPI.jl.
3. JuliaGeo: Focused on providing Julia interfaces to well-established external geospatial libraries. Key packages include GDAL.jl, Proj.jl, and LibGEOS.jl.

You can group geospatial data (geotables or domains) based on geometric criteria using the partition function. This function accepts a geotable or domain and a PartitionMethod object. The result is a lazy, indexable collection of the partitioned pieces, which can be accessed via indexing (e.g., Π[1]).

using GeoStats

# Create a sample geotable
gtb = georef((; z=rand(100, 100)))

# Partition using a specific method (e.g., BisectFractionPartition)
Π = partition(gtb, BisectFractionPartition((1.0, 1.0), fraction=0.7))

# Access the first partition
first_part = Π[1]

In the context of geospatial modeling, Gaussian process regression and Simple Kriging are essentially two names for the same concept.

However, there is a technical distinction: the derivation of Kriging estimators does not require distributional assumptions. For multivariate Gaussian distributions, Simple Kriging specifically provides the conditional expectation.

Geospatial data objects implement the Tables.jl interface, allowing them to behave like DataFrames.

• Attributes: Access columns using standard syntax (e.g., Ω.z or Ω[1, "column_name"]).
• Geometry: The geometry column is created on the fly to optimize performance. You can access the full table of attributes using values(Ω) and the underlying geospatial domain using domain(Ω).
• Rows: Access specific rows using Ω[row_index, :].

# Accessing data
Ω[1,:]          # First row
Ω[1,"geometry"] # Geometry of first row
Ω.z             # Column 'z'
values(Ω)        # Table of attributes
domain(Ω)        # Geospatial domain

Geostatistical interpolation models predict variables over geospatial domains. You can use these models within the Interpolate and InterpolateNeighbors transforms. These transforms support advanced options for change of support, probabilistic prediction, and neighborhood search.

All models operate on general Hilbert spaces, allowing you to interpolate any data type that implements scalar multiplication, vector addition, and an inner product.

using GeoStats

# Example setup
gtb = georef((; z=[1.,0.,1.]), [(25.,25.), (50.,75.), (75.,50.)])
grid = CartesianGrid(100, 100)

# Perform interpolation
itp = gtb |> Interpolate(grid, model=NN())

GeoStats.jl provides a hierarchy of transforms that interact with geospatial data in different ways. Understanding these categories is essential for choosing the right tool for your data:

1. Feature Transforms (TableTransforms.FeatureTransform): These operate exclusively on the attribute table (the data values) without changing the underlying geometry.
2. Geometric Transforms (Meshes.GeometricTransform): These operate on the geospatial domain (the geometry/coordinates) without changing the attribute values.
3. Geostatistical Transforms: These operate on both the geometries and the features simultaneously (e.g., Detrend, Potrace).

All transforms and pipelines implement a standard set of functions for application and reversal, including apply, revert, inverse, isrevertible, and isinvertible.

In GeoStats.jl, a "raster" is a Julia array that is georeferenced over a geospatial grid. This model is implemented via GeoTables over a Grid.

Key benefits include:

1. Convenient Syntax: Access grid elements using n-dimensional array syntax (e.g., slicing by indices or coordinate ranges).
2. Lazy Loading: Support for large datasets via packages like DiskArrays.jl, allowing you to work with data without loading the entire set into RAM.

GeoStats.jl integrates with the Makie.jl ecosystem for scientific visualization. To use these features, ensure a Makie backend is loaded in your session.

Key visualization functions include:

• viz: A recipe that takes a geospatial domain or a vector of geometries as input.
• viz!: An in-place version of the viz recipe.
• viewer: Takes a geotable as input and automatically calls viz on all columns that can be converted to colors by the Colorfy.jl package. It also automatically adds a color bar (cbar).
• cbar: Adds a color bar to a plot.

using GeoStats
import CairoMakie as Mke

# Example setup (assuming a Makie backend is loaded)
# viz(geodata)
# viewer(geotable)

GeoStats.jl defines two primary types of spatial statistical problems:

• Estimation: Given spatial data and a domain, the goal is to estimate variable(s) at unseen locations. This process typically provides a variance (or uncertainty) map alongside the estimates.
• Simulation: Given a domain and optionally spatial data, the goal is to generate multiple realizations of variable(s) that match the existing data where present.

GeoStats.jl provides two main approaches to clustering geospatial data into homogeneous regions:

1. Classical Clustering: These methods consider only the values of the geotable. They partition data based on feature similarity without explicitly accounting for the spatial domain (coordinates).
2. Geostatistical Clustering: These methods consider both the values and the domain of the geotable. They aim to find clusters that are both feature-homogeneous and spatially coherent.

To use these methods, you typically pipe a geotable into the desired clustering function.

using GeoStats

# Example geotable with values
gtb = georef((z=[10sin(i/10) + j for i in 1:4:100, j in 1:4:100],))

# Classical approach (values only)
ctb_classical = gtb |> KMedoids(5)

# Geostatistical approach (values + domain)
ctb_geostat = gtb |> GHC(5, 1.0)