If you find GeoStats.jl useful in your work, please consider citing it:
@ARTICLE{Hoffimann2018,
+ title={GeoStats.jl – High-performance geostatistics in Julia},
+ author={Hoffimann, Júlio},
+ journal={Journal of Open Source Software},
+ publisher={The Open Journal},
+ volume={3},
+ pages={692},
+ number={24},
+ ISSN={2475-9066},
+ DOI={10.21105/joss.00692},
+ url={https://dx.doi.org/10.21105/joss.00692},
+ year={2018},
+ month={Apr}
+}Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
We use the Zulip platform to chat and help the community of users. Consider creating an account and pressing ? on the keyboard there for navigation instructions.
Click on the image to join the channel:
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
The GeoStats.jl project is licensed under the MIT license:
MIT License
+
+Copyright (c) 2015 Júlio Hoffimann <julio.hoffimann@gmail.com> and contributors
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
+ ${display_result} +
+First off, thank you for considering contributing to GeoStats.jl. It’s people like you that make this project so much fun. Below are a few suggestions to speed up the collaboration process:
Contributing to an open-source project for the very first time can be a very daunting task. To make the process easier and more GitHub-beginner-friendly, the community has written an article about how to start contributing to open-source and overcome the mental and technical barriers that come associated with it. The article will also take you through the steps required to make your first contribution in detail.
If you are experiencing issues or have discovered a bug, please report it on GitHub. To make the resolution process easier, please include the version of Julia and GeoStats.jl in your writeup. These can be found with two commands:
julia> versioninfo()
+julia> using Pkg; Pkg.status()If you have suggestions of improvement or algorithms that you would like to see implemented, please open an issue on GitHub. Suggestions as well as feature requests are very welcome.
If you have code that you would like to contribute that is awesome! Please open an issue or reach out in our community channel before you create the pull request on GitHub so that we make sure your idea is aligned with the project goals.
After your idea is revised by project maintainers, you implement it online on Github or offline on your machine.
If the changes to the code are minimal, we recommend pressing . on the keyboard on any file in the GitHub repository of interest. This will open the VSCode editor on the web browser where you can implement the changes, commit them and submit a pull request.
If the changes require additional investigation and tests, please get the development version of the project by typing the following in the package manager:
] activate @geoThis will create a fresh environment called @geo where you can edit the project modules without effects on your global user environment. Next, go ahead and ask the package manager to develop the package of interest (e.g. GeoStatsFunctions.jl):
] dev GeoStatsFunctionsYou can modify the source code that was cloned in the .julia/dev folder and submit a pull request on GitHub later.
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
Given a table or array containing data, we can georeference these objects onto a geospatial domain with the georef function. The result of the georef function is a GeoTable. If you would like to learn this concept in more depth, check out the recording of our JuliaEO 2024 workshop:
+ +
GeoTables.georef — Functiongeoref(table, domain)Georeference table on domain from Meshes.jl
Examples
julia> georef((a=rand(100), b=rand(100)), CartesianGrid(10, 10))georef(table, geoms)Georeference table on vector of geometries geoms from Meshes.jl
Examples
julia> georef((a=rand(10), b=rand(10)), rand(Point, 10))georef(table, coords; [crs], [lenunit])Georeference table using coordinates coords of points.
Optionally, specify the coordinate reference system crs, which is set by default based on heuristics, and the lenunit (default to meters for unitless values) that can only be used in CRS types that allow this flexibility. Any CRS or EPSG/ESRI code from CoordRefSystems.jl is supported.
Examples
julia> georef((a=[1, 2, 3], b=[4, 5, 6], [(0, 0), (1, 1), (2, 2)])
+julia> georef((a=[1, 2, 3], b=[4, 5, 6], [(0, 0), (1, 1), (2, 2)], crs=LatLon)
+julia> georef((a=[1, 2, 3], b=[4, 5, 6], [(0, 0), (1, 1), (2, 2)], crs=EPSG{4326})
+julia> georef((a=[1, 2, 3], b=[4, 5, 6], [(0, 0), (1, 1), (2, 2)], lenunit=u"cm")georef(table, names; [crs], [lenunit])Georeference table using coordinates of points stored in column names.
Optionally, specify the coordinate reference system crs, which is set by default based on heuristics, and the lenunit (default to meters for unitless values) that can only be used in CRS types that allow this flexibility. Any CRS or EPSG/ESRI code from CoordRefSystems.jl is supported.
Examples
georef((a=rand(10), x=rand(10), y=rand(10)), ("x", "y"))
+georef((a=rand(10), x=rand(10), y=rand(10)), ("x", "y"), crs=LatLon)
+georef((a=rand(10), x=rand(10), y=rand(10)), ("x", "y"), crs=EPSG{4326})
+georef((a=rand(10), x=rand(10), y=rand(10)), ("x", "y"), lenunit=u"cm")georef(tuple)Georeference a named tuple on CartesianGrid(dims), with dims obtained from the arrays stored in the tuple.
Examples
julia> georef((a=rand(10, 10), b=rand(10, 10))) # 2D grid
+julia> georef((a=rand(10, 10, 10), b=rand(10, 10, 10))) # 3D gridThe functions values and domain can be used to retrieve the table of attributes and the underlying geospatial domain:
Base.values — Functionvalues(geotable, [rank])Return the values of geotable for a given rank as a table.
The rank is a non-negative integer that specifies the parametric dimension of the geometries of interest:
If the rank is not specified, it is assumed to be the rank of the elements of the domain.
GeoTables.domain — Functiondomain(geotable)Return underlying domain of the geotable.
using GeoStats
+import CairoMakie as Mke
+
+# helper function for plotting two
+# variables named T and P side by side
+function plot(data)
+ fig = Mke.Figure(size = (800, 400))
+ viz(fig[1,1], data.geometry, color = data.T)
+ viz(fig[1,2], data.geometry, color = data.P)
+ fig
+endplot (generic function with 1 method)Consider a table (e.g. DataFrame) with 25 samples of temperature T and pressure P:
using DataFrames
+
+table = DataFrame(T=rand(25), P=rand(25))| Row | T | P |
|---|---|---|
| Float64 | Float64 | |
| 1 | 0.477446 | 0.852434 |
| 2 | 0.944738 | 0.81652 |
| 3 | 0.704505 | 0.732916 |
| 4 | 0.936233 | 0.0320072 |
| 5 | 0.653428 | 0.393943 |
| 6 | 0.879168 | 0.694933 |
| 7 | 0.759295 | 0.601801 |
| 8 | 0.421052 | 0.155567 |
| 9 | 0.244843 | 0.782078 |
| 10 | 0.405972 | 0.767578 |
| 11 | 0.364723 | 0.750821 |
| 12 | 0.772122 | 0.362319 |
| 13 | 0.11324 | 0.625441 |
| 14 | 0.974632 | 0.373023 |
| 15 | 0.889477 | 0.15257 |
| 16 | 0.0433758 | 0.00187064 |
| 17 | 0.938691 | 0.205616 |
| 18 | 0.345166 | 0.800411 |
| 19 | 0.401956 | 0.286133 |
| 20 | 0.150285 | 0.880893 |
| 21 | 0.0752978 | 0.780687 |
| 22 | 0.23104 | 0.541633 |
| 23 | 0.545222 | 0.338354 |
| 24 | 0.663735 | 0.621553 |
| 25 | 0.698587 | 0.608446 |
We can georeference this table based on a given set of points:
georef(table, rand(Point, 25)) |> plot
or alternatively, georeference it on a 5x5 regular grid (5x5 = 25 samples):
georef(table, CartesianGrid(5, 5)) |> plot
Another common pattern in geospatial data is when the coordinates of the samples are already part of the table as columns. In this case, we can specify the column names:
table = DataFrame(T=rand(25), P=rand(25), X=rand(25), Y=rand(25), Z=rand(25))
+
+georef(table, ("X", "Y", "Z")) |> plot
Consider arrays (e.g. images) with data for various geospatial variables. We can georeference these arrays using a named tuple, and the framework will understand that the shape of the arrays should be preserved in a CartesianGrid:
T, P = rand(5, 5), rand(5, 5)
+
+georef((T=T, P=P)) |> plot
Alternatively, we can interpret the entries of the named tuple as columns in a table:
georef((T=vec(T), P=vec(P)), rand(Point, 25)) |> plot
The GeoIO.jl package can be used to load/save geospatial data from/to various file formats:
GeoIO.load — Functionload(fname, repair=true, layer=0, lenunit=m, kwargs...)Load geospatial table from file fname stored in any format.
Various repairs are performed on the stored geometries by default, including fixes of orientation in rings of polygons, removal of zero-area triangles, etc.
Some of the repairs can be expensive on large data sets. In that case, we recommend setting repair=false. Custom repairs can be performed with the Repair transform from Meshes.jl.
Optionally, specify the layer to read within the file, and the length unit lenunit of the coordinates when the format does not include units in its specification. Other kwargs are forwarded to the backend packages.
Please use the formats function to list all supported file formats.
Options
OFF
defaultcolor: default color of the geometries if the file does not have this data (default to RGBA(0.666, 0.666, 0.666, 0.666));CSV
coords: names of the columns with point coordinates (required option);CSV.File, see the CSV.jl documentation for more details;VTK formats (.vtu, .vtp, .vtr, .vts, .vti)
:MASK). If the column does not exist in the file, the full grid is returned;Common Data Model formats (NetCDF, GRIB)
x: name of the column with x coordinates (default to "x", "X", "lon", or "longitude");y: name of the column with y coordinates (default to "y", "Y", "lat", or "latitude");z: name of the column with z coordinates (default to "z", "Z", "depth", or "height");t: name of the column with time measurements (default to "t", "time", or "TIME");GeoJSON
numbertype: number type of geometry coordinates (default to Float64)GeoJSON.read, see the GeoJSON.jl documentation for more details;GSLIB
GslibIO.load, see the GslibIO.jl documentation for more details;Shapefile
Shapefile.read, see the Shapefile.jl documentation for more details;GeoParquet
GeoParquet.read, see the GeoParquet.jl documentation for more details;GeoTIFF, GeoPackage, KML
ArchGDAL.read, see the ArchGDAL.jl documentation for more details;Examples
# load coordinates of geojson file as Float64 (default)
+GeoIO.load("file.geojson")
+# load coordinates of geojson file as Float32
+GeoIO.load("file.geojson", numbertype=Float32)GeoIO.save — Functionsave(fname, geotable; kwargs...)Save geotable to file fname of given format based on the file extension.
Other kwargs are forwarded to the backend packages.
Please use the formats function to list all supported file formats.
Options
OFF
color: name of the column with geometry colors, if nothing the geometries will be saved without colors (default to nothing);MSH
vcolumn: name of the column in vertex table with node data, if nothing the geometries will be saved without node data (default to nothing);ecolumn: name of the column in element table with element data, if nothing the geometries will be saved without element data (default to nothing);STL
ascii: defines whether the file will be saved in ASCII format, otherwise Binary format will be used (default to false);CSV
coords: names of the columns where the point coordinates will be saved (default to "x", "y", "z");floatformat: C-style format string for float values (default to no formatting);CSV.write, see the CSV.jl documentation for more details;NetCDF
x: name of the column where the coordinate x will be saved (default to CRS coordinate name);y: name of the column where the coordinate y will be saved (default to CRS coordinate name);z: name of the column where the coordinate z will be saved (default to CRS coordinate name);t: name of the column where the time measurements will be saved (default to "t");GeoTIFF
options: list with options that will be passed to GDAL;GeoPackage
layername: name of the layer where the data will be saved (default to "data");options: dictionary with options that will be passed to GDAL;GSLIB
GslibIO.save, see the GslibIO.jl documentation for more details;Shapefile
Shapefile.write, see the Shapefile.jl documentation for more details;GeoJSON
GeoJSON.write, see the GeoJSON.jl documentation for more details;GeoParquet
GeoParquet.write, see the GeoParquet.jl documentation for more details;Examples
# overwrite an existing shapefile
+GeoIO.save("file.shp", force = true)GeoIO.formats — Functionformats([io]; sortby=:format)Displays in io (defaults to stdout if io is not given) a table with all formats supported by GeoIO.jl and the packages used to load and save each of them.
Optionally, sort the table by the :extension, :load or :save columns using the sortby argument.
using GeoIO
+
+zone = GeoIO.load("data/zone.shp")| PERIMETER | +ACRES | +MACROZONA | +Hectares | +area_m2 | +geometry | +
|---|---|---|---|---|---|
| Continuous | +Continuous | +Categorical | +Continuous | +Continuous | +MultiPolygon | +
| [NoUnits] | +[NoUnits] | +[NoUnits] | +[NoUnits] | +[NoUnits] | +🖈 GeodeticLatLon{SAD69} | +
| 5.8508e6 | +3.23145e7 | +Estuario | +1.30772e7 | +1.30772e11 | +Multi(21×PolyArea) | +
| 9.53947e6 | +2.50594e8 | +Fronteiras Antigas | +1.01412e8 | +1.01412e12 | +Multi(1×PolyArea) | +
| 1.01743e7 | +2.75528e8 | +Fronteiras Intermediarias | +1.11502e8 | +1.11502e12 | +Multi(1×PolyArea) | +
| 7.09612e6 | +1.61293e8 | +Fronteiras Novas | +6.5273e7 | +6.5273e11 | +Multi(2×PolyArea) | +
The GeoArtifacts.jl package provides utility functions to automatically download geospatial data from repositories on the internet.
using GeoArtifacts
+
+# download artifacts from naturalearthdata.com
+earth = NaturalEarth.naturalearth1("water")
+borders = NaturalEarth.borders()
+airports = NaturalEarth.airports()
+ports = NaturalEarth.ports()
+
+# initialize viewer with a coarse "raster"
+earth |> Upscale(10, 5) |> viewer
+
+# add other elements to the visualization
+viz!(borders.geometry, color = "cyan")
+viz!(airports.geometry, color = "black", pointsize=4, pointmarker='✈')
+viz!(ports.geometry, color = "blue", pointsize=4, pointmarker='⚓')
+
+# display current figure
+Mke.current_figure()
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
Declustered statistics are statistics that make use of geospatial coordinates in an attempt to correct potential sampling bias:
+ +
The following statistics have geospatial semantics:
Statistics.mean — Methodmean(data, v)
+mean(data, v, s)Declustered mean of geospatial data. Optionally, specify the variable v and the block side s.
Statistics.var — Methodvar(data, v)
+var(data, v, s)Declustered variance of geospatial data. Optionally, specify the variable v and the block side s.
Statistics.quantile — Methodquantile(data, v, p)
+quantile(data, v, p, s)Declustered quantile of geospatial data at probability p. Optionally, specify the variable v and the block side s.
A histogram with geospatial semantics is also available where the heights of the bins are adjusted based on the coordinates of the samples:
GeoStatsBase.EmpiricalHistogram — TypeEmpiricalHistogram(sdata, var, [s]; kwargs...)Spatial histogram of variable var in spatial data sdata. Optionally, specify the block side s and the keyword arguments kwargs for the fit(Histogram, ...) call.
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
A geospatial domain is a (discretized) region in physical space. For example, a collection of rain gauge stations can be represented as a PointSet in the map. Similarly, a collection of states in a given country can be represented as a GeometrySet of polygons.
We provide flexible domain types for advanced geospatial workflows via the Meshes.jl submodule. The domains can consist of sets (or "soups") of geometries as it is most common in GIS or be actual meshes with topological information.
Meshes.Domain — TypeDomain{M,CRS}A domain is an indexable collection of geometries (e.g. mesh) in a given manifold M with point coordinates specified in a coordinate reference system CRS.
Meshes.Mesh — TypeMesh{M,CRS,TP}A mesh of geometries in a given manifold M with point coordinates specified in a coordinate reference system CRS. Unlike a general domain, a mesh has a well-defined topology TP.
Meshes.PointSet — TypePointSet(points)A set of points (a.k.a. point cloud) representing a Domain.
Examples
All point sets below are the same and contain two points in R³:
julia> PointSet([Point(1,2,3), Point(4,5,6)])
+julia> PointSet(Point(1,2,3), Point(4,5,6))
+julia> PointSet([(1,2,3), (4,5,6)])
+julia> PointSet((1,2,3), (4,5,6))pset = PointSet(rand(Point, 100))
+
+viz(pset)
Meshes.GeometrySet — TypeGeometrySet(geometries)A set of geometries representing a Domain.
Examples
Set containing two balls centered at (0.0, 0.0) and (1.0, 1.0):
julia> GeometrySet([Ball((0.0, 0.0)), Ball((1.0, 1.0))])Notes
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)
Meshes.RegularGrid — TypeRegularGrid(dims, origin, spacing)A regular grid with dimensions dims, lower left corner at origin and cell spacing spacing. The three arguments must have the same length.
RegularGrid(dims, origin, spacing, offset)A regular grid with dimensions dims, with lower left corner of element offset at origin and cell spacing spacing.
RegularGrid(start, finish, dims=dims)Alternatively, construct a regular grid from a start point to a finish with dimensions dims.
RegularGrid(start, finish, spacing)Alternatively, construct a regular grid from a start point to a finish point using a given spacing.
Examples
RegularGrid((10, 20), Point(LatLon(30.0°, 60.0°)), (1.0, 1.0)) # add coordinate units to spacing
+RegularGrid((10, 20), Point(Polar(0.0cm, 0.0rad)), (10.0mm, 1.0rad)) # convert spacing units to coordinate units
+RegularGrid((10, 20), Point(Mercator(0.0, 0.0)), (1.5, 1.5))
+RegularGrid((10, 20, 30), Point(Cylindrical(0.0, 0.0, 0.0)), (3.0, 2.0, 1.0))See also CartesianGrid.
Meshes.CartesianGrid — TypeCartesianGrid(dims, origin, spacing)A Cartesian grid with dimensions dims, lower left corner at origin and cell spacing spacing. The three arguments must have the same length.
CartesianGrid(dims, origin, spacing, offset)A Cartesian grid with dimensions dims, with lower left corner of element offset at origin and cell spacing spacing.
CartesianGrid(start, finish, dims=dims)Alternatively, construct a Cartesian grid from a start point (lower left) to a finish point (upper right).
CartesianGrid(start, finish, spacing)Alternatively, construct a Cartesian grid from a start point to a finish point using a given spacing.
CartesianGrid(dims)
+CartesianGrid(dim1, dim2, ...)Finally, a Cartesian grid can be constructed by only passing the dimensions dims as a tuple, or by passing each dimension dim1, dim2, ... separately. In this case, the origin and spacing default to (0,0,...) and (1,1,...).
CartesianGrid is an alias to RegularGrid with Cartesian CRS.
Examples
Create a 3D grid with 100x100x50 hexahedrons:
julia> CartesianGrid(100, 100, 50)Create a 2D grid with 100 x 100 quadrangles and origin at (10.0, 20.0):
julia> CartesianGrid((100, 100), (10.0, 20.0), (1.0, 1.0))Create a 1D grid from -1 to 1 with 100 segments:
julia> CartesianGrid((-1.0,), (1.0,), dims=(100,))See also RegularGrid.
grid = CartesianGrid(10, 10, 10)
+
+viz(grid, showsegments = true)
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
GeoStats — Module
GeoStats.jl is an extensible framework for geospatial data science and geostatistical modeling fully written in Julia. It is comprised of several modules for advanced geometric processing, state-of-the-art geostatistical algorithms and sophisticated visualization of geospatial data.
All further information is provided in the online documentation.
If you have found this software useful, please consider starring it on GitHub. This gives us an accurate lower bound of the (satisfied) user count.
Organizations using the framework:
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+
+
+ 
+
+
Would like to become a sponsor? Press the sponsor button in our GitHub repository.
In many fields of science, such as mining engineering, hydrogeology, petroleum engineering, and environmental sciences, traditional statistical methods fail to provide unbiased estimates of resources due to the presence of geospatial correlation. Geostatistics (a.k.a. geospatial statistics) is the branch of statistics developed to overcome this limitation. Particularly, it is the branch that takes geospatial coordinates of data into account. Some major highlights of GeoStats.jl are:
Our JuliaCon2021 talk provides an overview of geostatistical learning, which is one of the many geostatistical problems addressed by the software:
+ +
Consider reading the Geospatial Data Science with Julia book before reading this documentation. If you have questions, or would like to brainstorm ideas, don't hesitate to start a topic in our community channel.
+
+ 
+
+
Get the latest stable release with Julia's package manager:
using Pkg
+Pkg.add("GeoStats")Below is an example of geostatistical interpolation of point data over a Cartesian grid with a Kriging model:
# 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)
+
+# visualize the solution
+viz(interp.geometry, color = interp.Z)
The project is split into various packages:
| Package | Description |
|---|---|
| GeoStats.jl | Main package reexporting full stack of packages for geostatistics. |
| Meshes.jl | Computational geometry and advanced meshing algorithms. |
| CoordRefSystems.jl | Unitful coordinate reference systems. |
| GeoTables.jl | Geospatial tables compatible with the framework. |
| DataScienceTraits.jl | Traits for geospatial data science. |
| TableTransforms.jl | Transforms and pipelines with tabular data. |
| StatsLearnModels.jl | Statistical learning models for geospatial prediction. |
| GeoStatsBase.jl | Base package with core geostatistical definitions. |
| GeoStatsFunctions.jl | Geostatistical functions and related tools. |
| GeoStatsModels.jl | Geostatistical models for geospatial interpolation. |
| GeoStatsProcesses.jl | Geostatistical processes for geospatial simulation. |
| GeoStatsTransforms.jl | Geostatistical transforms for geospatial data. |
| GeoStatsValidation.jl | Geostatistical validation methods. |
Other packages can be installed separately for additional functionality:
| Package | Description |
|---|---|
| GeoIO.jl | Load/save geospatial tables in various formats. |
| GeoArtifacts.jl | Artifacts for geospatial data science. |
| DrillHoles.jl | Desurvey/composite drillhole data. |
If you find this software useful in your work, please consider citing it:
@ARTICLE{Hoffimann2018,
+ title={GeoStats.jl – High-performance geostatistics in Julia},
+ author={Hoffimann, Júlio},
+ journal={Journal of Open Source Software},
+ publisher={The Open Journal},
+ volume={3},
+ pages={692},
+ number={24},
+ ISSN={2475-9066},
+ DOI={10.21105/joss.00692},
+ url={https://dx.doi.org/10.21105/joss.00692},
+ year={2018},
+ month={Apr}
+}We ❤ to see our list of publications growing.
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
This section describes the Kriging models used in the Interpolate and InterpolateNeighbors transforms, which provide options for neighborhood search, change of support, etc.
A Kriging model has the form:
\[\hat{Z}(\x_0) = \lambda_1 Z(\x_1) + \lambda_2 Z(\x_2) + \cdots + \lambda_n Z(\x_n),\quad \x_i \in \R^m, \lambda_i \in \R\]
with $Z\colon \R^m \times \Omega \to \R$ a random field.
This package implements the following Kriging variants:
All these variants follow the same interface: an object is first created with a given set of parameters, it is then combined with the data to obtain predictions at new geometries.
The fit function takes care of building the Kriging system and factorizing the LHS with an appropriate decomposition (e.g. Cholesky, LU), and the predict (or predictprob) function performs the estimation for a given variable and geometry.
All variants work with general Hilbert spaces, meaning that one can interpolate any data type that implements scalar multiplication, vector addition and inner product.
In Simple Kriging, the mean $\mu$ of the random field is assumed to be constant and known. The resulting linear system is:
\[\begin{bmatrix} +cov(\x_1,\x_1) & cov(\x_1,\x_2) & \cdots & cov(\x_1,\x_n) \\ +cov(\x_2,\x_1) & cov(\x_2,\x_2) & \cdots & cov(\x_2,\x_n) \\ +\vdots & \vdots & \ddots & \vdots \\ +cov(\x_n,\x_1) & cov(\x_n,\x_2) & \cdots & cov(\x_n,\x_n) +\end{bmatrix} +\begin{bmatrix} +\lambda_1 \\ +\lambda_2 \\ +\vdots \\ +\lambda_n +\end{bmatrix} += +\begin{bmatrix} +cov(\x_1,\x_0) \\ +cov(\x_2,\x_0) \\ +\vdots \\ +cov(\x_n,\x_0) +\end{bmatrix}\]
or in matricial form $\C\l = \c$. We subtract the given mean from the observations $\boldsymbol{y} = \z - \mu \1$ and compute the mean and variance at location $\x_0$:
\[\mu(\x_0) = \mu + \boldsymbol{y}^\top \l\]
\[\sigma^2(\x_0) = cov(0) - \c^\top \l\]
GeoStatsModels.SimpleKriging — TypeSimpleKriging(γ, μ)Simple Kriging with variogram model γ and constant mean μ.
Notes
In Ordinary Kriging the mean of the random field is assumed to be constant and unknown. The resulting linear system is:
\[\begin{bmatrix} +\G & \1 \\ +\1^\top & 0 +\end{bmatrix} +\begin{bmatrix} +\l \\ +\nu +\end{bmatrix} += +\begin{bmatrix} +\g \\ +1 +\end{bmatrix}\]
with $\nu$ the Lagrange multiplier associated with the constraint $\1^\top \l = 1$. The mean and variance at location $\x_0$ are given by:
\[\mu(\x_0) = \z^\top \l\]
\[\sigma^2(\x_0) = \begin{bmatrix} \g \\ 1 \end{bmatrix}^\top \begin{bmatrix} \l \\ \nu \end{bmatrix}\]
GeoStatsModels.OrdinaryKriging — TypeOrdinaryKriging(γ)Ordinary Kriging with variogram model γ.
In Universal Kriging, the mean of the random field is assumed to be a polynomial of the geospatial coordinates:
\[\mu(\x) = \sum_{k=1}^{N_d} \beta_k f_k(\x)\]
with $N_d$ monomials $f_k$ of degree up to $d$. For example, in 2D there are $6$ monomials of degree up to $2$:
\[\mu(x_1,x_2) = \beta_1 1 + \beta_2 x_1 + \beta_3 x_2 + \beta_4 x_1 x_2 + \beta_5 x_1^2 + \beta_6 x_2^2\]
The choice of the degree $d$ determines the size of the polynomial matrix
\[\F = +\begin{bmatrix} +f_1(\x_1) & f_2(\x_1) & \cdots & f_{N_d}(\x_1) \\ +f_1(\x_2) & f_2(\x_2) & \cdots & f_{N_d}(\x_2) \\ +\vdots & \vdots & \ddots & \vdots \\ +f_1(\x_n) & f_2(\x_n) & \cdots & f_{N_d}(\x_n) +\end{bmatrix}\]
and polynomial vector $\f = \begin{bmatrix} f_1(\x_0) & f_2(\x_0) & \cdots & f_{N_d}(\x_0) \end{bmatrix}^\top$.
The variogram determines the variogram matrix:
\[\G = +\begin{bmatrix} +\gamma(\x_1,\x_1) & \gamma(\x_1,\x_2) & \cdots & \gamma(\x_1,\x_n) \\ +\gamma(\x_2,\x_1) & \gamma(\x_2,\x_2) & \cdots & \gamma(\x_2,\x_n) \\ +\vdots & \vdots & \ddots & \vdots \\ +\gamma(\x_n,\x_1) & \gamma(\x_n,\x_2) & \cdots & \gamma(\x_n,\x_n) +\end{bmatrix}\]
and the variogram vector $\g = \begin{bmatrix} \gamma(\x_1,\x_0) & \gamma(\x_2,\x_0) & \cdots & \gamma(\x_n,\x_0) \end{bmatrix}^\top$.
The resulting linear system is:
\[\begin{bmatrix} +\G & \F \\ +\F^\top & \boldsymbol{0} +\end{bmatrix} +\begin{bmatrix} +\l \\ +\boldsymbol{\nu} +\end{bmatrix} += +\begin{bmatrix} +\g \\ +\f +\end{bmatrix}\]
with $\boldsymbol{\nu}$ the Lagrange multipliers associated with the universal constraints. The mean and variance at location $\x_0$ are given by:
\[\mu(\x_0) = \z^\top \l\]
\[\sigma^2(\x_0) = \begin{bmatrix}\g \\ \f\end{bmatrix}^\top \begin{bmatrix}\l \\ \boldsymbol{\nu}\end{bmatrix}\]
GeoStatsModels.UniversalKriging — TypeUniversalKriging(γ, degree, dim)Universal Kriging with variogram model γ and polynomial of given degree on dim coordinates.
Notes
OrdinaryKriging is recovered for 0th degree polynomialExternalDriftKrigingIn External Drift Kriging, the mean of the random field is assumed to be a combination of known smooth functions:
\[\mu(\x) = \sum_k \beta_k m_k(\x)\]
Differently than Universal Kriging, the functions $m_k$ are not necessarily polynomials of the geospatial coordinates. In practice, they represent a list of variables that is strongly correlated (and co-located) with the variable being estimated.
External drifts are known to cause numerical instability. Give preference to other Kriging variants if possible.
GeoStatsModels.ExternalDriftKriging — TypeExternalDriftKriging(γ, drifts)External Drift Kriging with variogram model γ and external drifts functions. A drift function p -> v maps a point p to a value v.
Notes
p -> 1) for unbiased estimationOrdinaryKriging is recovered for drifts = [p -> 1]UniversalKrigingSettings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
Below is the list of types and functions mentioned in the documentation.
GeoStatsBase.DatamineAnglesGeoStatsBase.EmpiricalHistogramGeoStatsBase.GslibAnglesGeoStatsBase.MinesightAnglesGeoStatsBase.VulcanAnglesGeoStatsFunctions.CircularVariogramGeoStatsFunctions.CubicVariogramGeoStatsFunctions.EmpiricalTransiogramGeoStatsFunctions.EmpiricalTransioplaneGeoStatsFunctions.EmpiricalVariogramGeoStatsFunctions.EmpiricalVarioplaneGeoStatsFunctions.ExponentialTransiogramGeoStatsFunctions.ExponentialVariogramGeoStatsFunctions.GaussianTransiogramGeoStatsFunctions.GaussianVariogramGeoStatsFunctions.LinearTransiogramGeoStatsFunctions.MaternVariogramGeoStatsFunctions.MatrixExponentialTransiogramGeoStatsFunctions.NuggetEffectGeoStatsFunctions.PentasphericalVariogramGeoStatsFunctions.PiecewiseLinearTransiogramGeoStatsFunctions.PowerVariogramGeoStatsFunctions.SineHoleVariogramGeoStatsFunctions.SphericalTransiogramGeoStatsFunctions.SphericalVariogramGeoStatsFunctions.WeightedLeastSquaresGeoStatsModels.ExternalDriftKrigingGeoStatsModels.OrdinaryKrigingGeoStatsModels.SimpleKrigingGeoStatsModels.UniversalKrigingGeoStatsProcesses.BinomialProcessGeoStatsProcesses.ClusterProcessGeoStatsProcesses.EnsembleGeoStatsProcesses.FFTMethodGeoStatsProcesses.GaussianProcessGeoStatsProcesses.InhibitionProcessGeoStatsProcesses.LUMethodGeoStatsProcesses.LindgrenProcessGeoStatsProcesses.PoissonProcessGeoStatsProcesses.QuiltingProcessGeoStatsProcesses.RandomThinningGeoStatsProcesses.SEQMethodGeoStatsProcesses.StrataProcessGeoStatsProcesses.TuringProcessGeoStatsTransforms.CookieCutterGeoStatsTransforms.DetrendGeoStatsTransforms.GHCGeoStatsTransforms.GSCGeoStatsTransforms.InterpolateGeoStatsTransforms.InterpolateMissingGeoStatsTransforms.InterpolateNaNGeoStatsTransforms.InterpolateNeighborsGeoStatsTransforms.PotraceGeoStatsTransforms.RasterizeGeoStatsTransforms.SLICGeoStatsTransforms.SimulateGeoStatsTransforms.UniqueCoordsGeoStatsValidation.BlockValidationGeoStatsValidation.DensityRatioValidationGeoStatsValidation.KFoldValidationGeoStatsValidation.LeaveBallOutGeoStatsValidation.LeaveOneOutGeoStatsValidation.WeightedValidationMeshes.CartesianGridMeshes.DomainMeshes.GeometricTransformMeshes.GeometrySetMeshes.MeshMeshes.PointSetMeshes.RegularGridStatsLearnModels.LearnTableTransforms.FeatureTransformTableTransforms.SelectBase.randBase.randBase.unionBase.valuesGeoIO.formatsGeoIO.loadGeoIO.saveGeoStatsFunctions.DirectionalTransiogramGeoStatsFunctions.DirectionalVariogramGeoStatsFunctions.PlanarTransiogramGeoStatsFunctions.PlanarVariogramGeoStatsFunctions.fitGeoStatsFunctions.isstationaryGeoStatsFunctions.structuresGeoStatsFunctions.transioplotGeoStatsFunctions.varioplotGeoStatsProcesses.ishomogeneousGeoStatsProcesses.thinGeoStatsValidation.cverrorGeoTables.domainGeoTables.geojoinGeoTables.georefGeoTables.tablejoinGeoTables.viewerMeshes.isisotropicMeshes.vizMeshes.viz!Statistics.meanStatistics.quantileStatistics.varTransformsBase.applyTransformsBase.inverseTransformsBase.isinvertibleTransformsBase.isrevertibleTransformsBase.reapplyTransformsBase.revertSettings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
We provide a geospatial version of the split-apply-combine pattern:
+ +
GeoTables.@groupby — Macro@groupby(geotable, col₁, col₂, ..., colₙ)
+@groupby(geotable, [col₁, col₂, ..., colₙ])
+@groupby(geotable, (col₁, col₂, ..., colₙ))Group geospatial geotable by columns col₁, col₂, ..., colₙ.
@groupby(geotable, regex)Group geospatial geotable by columns that match with regex.
Examples
@groupby(geotable, 1, 3, 5)
+@groupby(geotable, [:a, :c, :e])
+@groupby(geotable, ("a", "c", "e"))
+@groupby(geotable, r"[ace]")GeoTables.@transform — Macro@transform(geotable, :col₁ = expr₁, :col₂ = expr₂, ..., :colₙ = exprₙ)Returns geospatial geotable with columns col₁, col₂, ..., colₙ computed with expressions expr₁, expr₂, ..., exprₙ.
See also: @groupby.
Examples
@transform(geotable, :z = :x + 2*:y)
+@transform(geotable, :w = :x^2 - :y^2)
+@transform(geotable, :sinx = sin(:x), :cosy = cos(:y))
+
+groups = @groupby(geotable, :y)
+@transform(groups, :logx = log(:x))
+@transform(groups, :expz = exp(:z))
+
+@transform(geotable, {"z"} = {"x"} - 2*{"y"})
+xnm, ynm, znm = :x, :y, :z
+@transform(geotable, {znm} = {xnm} - 2*{ynm})GeoTables.@combine — Macro@combine(geotable, :col₁ = expr₁, :col₂ = expr₂, ..., :colₙ = exprₙ)Returns geospatial geotable with columns :col₁, :col₂, ..., :colₙ computed with reduction expressions expr₁, expr₂, ..., exprₙ.
If a reduction expression is not defined for the :geometry column, the geometries will be reduced using Multi.
See also: @groupby.
Examples
@combine(geotable, :x_sum = sum(:x))
+@combine(geotable, :x_mean = mean(:x))
+@combine(geotable, :x_mean = mean(:x), :geometry = centroid(:geometry))
+
+groups = @groupby(geotable, :y)
+@combine(groups, :x_prod = prod(:x))
+@combine(groups, :x_median = median(:x))
+@combine(groups, :x_median = median(:x), :geometry = centroid(:geometry))
+
+@combine(geotable, {"z"} = sum({"x"}) + prod({"y"}))
+xnm, ynm, znm = :x, :y, :z
+@combine(geotable, {znm} = sum({xnm}) + prod({ynm}))GeoTables.geojoin — Functiongeojoin(geotable₁, geotable₂, var₁ => agg₁, ..., varₙ => aggₙ; kind=:left, pred=intersects, on=nothing)Joins geotable₁ with geotable₂ using a certain kind of join and predicate function pred that takes two geometries and returns a boolean ((g1, g2) -> g1 ⊆ g2).
Optionally, add a variable value match to join, in addition to geometric match, by passing a single name or list of variable names (strings or symbols) to the on keyword argument. The variable value match use the isequal function.
Whenever two or more matches are encountered, aggregate varᵢ with aggregation function aggᵢ. If no aggregation function is provided for a variable, then the aggregation function will be selected according to the scientific types: mean for continuous and first otherwise.
Kinds
:left - Returns all rows of geotable₁ filling entries with missing when there is no match in geotable₂.:inner - Returns the subset of rows of geotable₁ that has a match in geotable₂.Examples
geojoin(gtb1, gtb2)
+geojoin(gtb1, gtb2, 1 => mean)
+geojoin(gtb1, gtb2, :a => mean, :b => std)
+geojoin(gtb1, gtb2, "a" => mean, pred=issubset)
+geojoin(gtb1, gtb2, on=:a)
+geojoin(gtb1, gtb2, kind=:inner, on=["a", "b"])See also tablejoin.
GeoTables.tablejoin — Functiontablejoin(geotable, table, var₁ => agg₁, ..., varₙ => aggₙ; kind=:left, on)Joins geotable with table using the values of the on variables to match rows with a certain kind of join.
on variables can be a single name or list of variable names (strings or symbols). The variable value match use the isequal function.
Whenever two or more matches are encountered, aggregate varᵢ with aggregation function aggᵢ. If no aggregation function is provided for a variable, then the aggregation function will be selected according to the scientific types: mean for continuous and first otherwise.
Kinds
:left - Returns all rows of geotable filling entries with missing when there is no match in table.:inner - Returns the subset of rows of geotable that has a match in table.Examples
tablejoin(gtb, tab, on=:a)
+tablejoin(gtb, tab, 1 => mean, on=:a)
+tablejoin(gtb, tab, :a => mean, :b => std, on=:a)
+tablejoin(gtb, tab, "a" => mean, on=[:a, :b])
+tablejoin(gtb, tab, kind=:inner, on=["a", "b"])See also geojoin.
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
A geostatistical workflow often consists of four steps:
In this section, we walk through these steps to illustrate some of the features of the project. In the case of geostatistical modeling, we will specifically explore geostatistical learning models. If you prefer learning from video, check out the recording of our JuliaEO 2023 workshop:
+ +
We assume that the following packages are loaded throughout the code examples:
using GeoStats
+import CairoMakie as MkeThe GeoStats.jl package exports the full stack for geospatial data science and geostatistical modeling. The CairoMakie.jl package is one of the possible visualization backends from the Makie.jl project.
If you are new to Julia and have never heard of Makie.jl before, here are a few tips to help you choose between the different backends:
WGLMakie.jl is the preferred backend for interactive visualization on the web browser. It integrates well with Pluto.jl notebooks and other web-based applications.
GLMakie.jl is the preferred backend for interactive high-performance visualization. It leverages native graphical resources and doesn't require a web browser to function.
CairoMakie.jl is the preferred backend for publication-quality static visualization. It requires less computing power and is therefore recommended for those users with modest laptops.
We recommend importing the Makie.jl backend as Mke to avoid polluting the Julia session with names from the visualization stack.
The Julia ecosystem for loading geospatial data is comprised of several low-level packages such as Shapefile.jl and GeoJSON.jl, which define their own very basic geometry types. Instead of requesting users to learn the so called GeoInterface.jl to handle these types, we provide the high-level GeoIO.jl package to load any file with geospatial data into well-tested geometries from the Meshes.jl submodule:
using GeoIO
+
+zone = GeoIO.load("data/zone.shp")| PERIMETER | +ACRES | +MACROZONA | +Hectares | +area_m2 | +geometry | +
|---|---|---|---|---|---|
| Continuous | +Continuous | +Categorical | +Continuous | +Continuous | +MultiPolygon | +
| [NoUnits] | +[NoUnits] | +[NoUnits] | +[NoUnits] | +[NoUnits] | +🖈 GeodeticLatLon{SAD69} | +
| 5.8508e6 | +3.23145e7 | +Estuario | +1.30772e7 | +1.30772e11 | +Multi(21×PolyArea) | +
| 9.53947e6 | +2.50594e8 | +Fronteiras Antigas | +1.01412e8 | +1.01412e12 | +Multi(1×PolyArea) | +
| 1.01743e7 | +2.75528e8 | +Fronteiras Intermediarias | +1.11502e8 | +1.11502e12 | +Multi(1×PolyArea) | +
| 7.09612e6 | +1.61293e8 | +Fronteiras Novas | +6.5273e7 | +6.5273e11 | +Multi(2×PolyArea) | +
Various functions are defined over these geometries, for instance:
sum(area, zone.geometry)2.9037868929915425e12 m^2sum(perimeter, zone.geometry)3.2680213406922124e7 mPlease check the Meshes.jl documentation for more details.
We highly recommend using Meshes.jl geometries in geospatial workflows as they were carefully designed to accomodate advanced features of the GeoStats.jl framework. Any other geometry type will likely fail with our geostatistical algorithms and pipelines.
Geospatial data can be derived from other Julia variables. For example, given a Julia array, which is not attached to any coordinate system, we can georeference the array using the georef function:
Z = [10sin(i/10) + 2j for i in 1:50, j in 1:50]
+
+Ω = georef((Z=Z,))| Z | +geometry | +
|---|---|
| Continuous | +Quadrangle | +
| [NoUnits] | +🖈 Cartesian{NoDatum} | +
| 2.99833 | +Quadrangle((x: 0.0 m, y: 0.0 m), ..., (x: 0.0 m, y: 1.0 m)) | +
| 3.98669 | +Quadrangle((x: 1.0 m, y: 0.0 m), ..., (x: 1.0 m, y: 1.0 m)) | +
| 4.9552 | +Quadrangle((x: 2.0 m, y: 0.0 m), ..., (x: 2.0 m, y: 1.0 m)) | +
| 5.89418 | +Quadrangle((x: 3.0 m, y: 0.0 m), ..., (x: 3.0 m, y: 1.0 m)) | +
| 6.79426 | +Quadrangle((x: 4.0 m, y: 0.0 m), ..., (x: 4.0 m, y: 1.0 m)) | +
| 7.64642 | +Quadrangle((x: 5.0 m, y: 0.0 m), ..., (x: 5.0 m, y: 1.0 m)) | +
| 8.44218 | +Quadrangle((x: 6.0 m, y: 0.0 m), ..., (x: 6.0 m, y: 1.0 m)) | +
| 9.17356 | +Quadrangle((x: 7.0 m, y: 0.0 m), ..., (x: 7.0 m, y: 1.0 m)) | +
| 9.83327 | +Quadrangle((x: 8.0 m, y: 0.0 m), ..., (x: 8.0 m, y: 1.0 m)) | +
| 10.4147 | +Quadrangle((x: 9.0 m, y: 0.0 m), ..., (x: 9.0 m, y: 1.0 m)) | +
| ⋮ | +⋮ | +
Default coordinates are assigned based on the size of the array, and different configurations can be obtained with different methods (see Data).
Geospatial data can be visualized with the viz recipe function:
viz(Ω.geometry, color = Ω.Z)
Alternatively, we provide a basic scientific viewer to visualize all viewable variables in the data with a colorbar and other interactive elements (see Visualization):
viewer(Ω)
Julia supports unicode symbols with $\LaTeX$ syntax. We can type \Omega and press TAB to get the symbol Ω in the example above. This autocompletion works in various text editors, including the VSCode editor with the Julia extension.
Our geospatial data implements the Tables.jl interface, which means that they can be accessed as if they were tables with samples in the rows and variables in the columns. In this case, a special column named geometry is created on the fly, row by row, containing Meshes.jl geometries.
For those familiar with the productive DataFrames.jl interface, there is nothing new:
Ω[1,:](Z = 2.9983341664682817, geometry = Quadrangle((x: 0.0 m, y: 0.0 m), ..., (x: 0.0 m, y: 1.0 m)))Ω[1,:geometry]Quadrangle
+├─ Point(x: 0.0 m, y: 0.0 m)
+├─ Point(x: 1.0 m, y: 0.0 m)
+├─ Point(x: 1.0 m, y: 1.0 m)
+└─ Point(x: 0.0 m, y: 1.0 m)Ω.Z2500-element Vector{Float64}:
+ 2.9983341664682817
+ 3.9866933079506124
+ 4.955202066613396
+ 5.894183423086505
+ 6.79425538604203
+ 7.6464247339503535
+ 8.44217687237691
+ 9.173560908995228
+ 9.833269096274833
+ 10.414709848078965
+ ⋮
+ 91.28424227586412
+ 90.83834063250545
+ 90.48397926110484
+ 90.22469882334903
+ 90.06308996366536
+ 90.000767424359
+ 90.0383539116416
+ 90.17547387375667
+ 90.41075725336862However, notice that our implementation performs some clever optimizations behind the scenes to avoid expensive creation of geometries:
Ω.geometry50×50 CartesianGrid
+├─ minimum: Point(x: 0.0 m, y: 0.0 m)
+├─ maximum: Point(x: 50.0 m, y: 50.0 m)
+└─ spacing: (1.0 m, 1.0 m)We can always retrieve the table of attributes (or features) with the function values and the underlying geospatial domain with the function domain. This can be useful for writing algorithms that are type-stable and depend purely on the feature values:
values(Ω)(Z = [2.9983341664682817, 3.9866933079506124, 4.955202066613396, 5.894183423086505, 6.79425538604203, 7.6464247339503535, 8.44217687237691, 9.173560908995228, 9.833269096274833, 10.414709848078965 … 91.8172288893559, 91.28424227586412, 90.83834063250545, 90.48397926110484, 90.22469882334903, 90.06308996366536, 90.000767424359, 90.0383539116416, 90.17547387375667, 90.41075725336862],)or on the geometries:
domain(Ω)50×50 CartesianGrid
+├─ minimum: Point(x: 0.0 m, y: 0.0 m)
+├─ maximum: Point(x: 50.0 m, y: 50.0 m)
+└─ spacing: (1.0 m, 1.0 m)It is easy to design advanced geospatial pipelines that operate on both the table of features and the underlying geospatial domain:
# 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
Coordinates before pipeline:
boundingbox(Ω.geometry)Box
+├─ min: Point(x: 0.0 m, y: 0.0 m)
+└─ max: Point(x: 50.0 m, y: 50.0 m)and after pipeline:
boundingbox(Ω̂.geometry)Box
+├─ min: Point(x: -0.5 m, y: -0.5 m)
+└─ max: Point(x: 0.5 m, y: 0.5 m)These pipelines are revertible meaning that we can transform the data, perform geostatistical modeling, and revert the pipelines to obtain estimates in the original sample space (see Transforms).
We provide three macros @groupby, @transform and @combine for powerful geospatial split-apply-combine patterns, as well as the function geojoin for advanced geospatial joins.
@transform(Ω, :W = 2 * :Z * area(:geometry))| Z | +W | +geometry | +
|---|---|---|
| Continuous | +Continuous | +Quadrangle | +
| [NoUnits] | +[m^2] | +🖈 Cartesian{NoDatum} | +
| 2.99833 | +5.99667 m^2 | +Quadrangle((x: 0.0 m, y: 0.0 m), ..., (x: 0.0 m, y: 1.0 m)) | +
| 3.98669 | +7.97339 m^2 | +Quadrangle((x: 1.0 m, y: 0.0 m), ..., (x: 1.0 m, y: 1.0 m)) | +
| 4.9552 | +9.9104 m^2 | +Quadrangle((x: 2.0 m, y: 0.0 m), ..., (x: 2.0 m, y: 1.0 m)) | +
| 5.89418 | +11.7884 m^2 | +Quadrangle((x: 3.0 m, y: 0.0 m), ..., (x: 3.0 m, y: 1.0 m)) | +
| 6.79426 | +13.5885 m^2 | +Quadrangle((x: 4.0 m, y: 0.0 m), ..., (x: 4.0 m, y: 1.0 m)) | +
| 7.64642 | +15.2928 m^2 | +Quadrangle((x: 5.0 m, y: 0.0 m), ..., (x: 5.0 m, y: 1.0 m)) | +
| 8.44218 | +16.8844 m^2 | +Quadrangle((x: 6.0 m, y: 0.0 m), ..., (x: 6.0 m, y: 1.0 m)) | +
| 9.17356 | +18.3471 m^2 | +Quadrangle((x: 7.0 m, y: 0.0 m), ..., (x: 7.0 m, y: 1.0 m)) | +
| 9.83327 | +19.6665 m^2 | +Quadrangle((x: 8.0 m, y: 0.0 m), ..., (x: 8.0 m, y: 1.0 m)) | +
| 10.4147 | +20.8294 m^2 | +Quadrangle((x: 9.0 m, y: 0.0 m), ..., (x: 9.0 m, y: 1.0 m)) | +
| ⋮ | +⋮ | +⋮ | +
These are very useful for geospatial data science as they hide the complexity of the geometry column. For more information, check the Queries section of the documentation.
Having defined the geospatial data objects, we proceed and define the geostatistical learning model. Let's assume that we have geopatial data with some variable that we want to predict in a supervised learning setting. We load the data from a CSV file, and inspect the available columns:
using CSV
+using DataFrames
+
+table = CSV.File("data/agriculture.csv") |> DataFrame
+
+first(table, 5)| Row | band1 | band2 | band3 | band4 | crop | polygon | x | y |
|---|---|---|---|---|---|---|---|---|
| Int64 | Int64 | Int64 | Int64 | Int64 | Float64 | Float64 | Float64 | |
| 1 | 121 | 91 | 70 | 384 | 1 | 222.0 | 746.0 | 717.0 |
| 2 | 112 | 82 | 58 | 336 | 2 | 264.0 | 953.0 | 291.0 |
| 3 | 120 | 86 | 55 | 389 | 1 | 363.0 | 1050.0 | 796.0 |
| 4 | 125 | 92 | 78 | 255 | 2 | 295.0 | 787.0 | 921.0 |
| 5 | 117 | 77 | 59 | 308 | 1 | 347.0 | 947.0 | 758.0 |
Columns band1, band2, ..., band4 represent four satellite bands for different locations (x, y) in this region. The column crop has the crop type for each location that was labeled manually with the purpose of fitting a learning model.
We can now georeference the table and plot some of the variables:
Ω = 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
Similar to a generic statistical learning workflow, we split the data into "train" and "test" sets. The main difference here is that our geospatial geosplit function accepts a separating plane specified by its normal direction (1,-1):
Ωs, Ωt = geosplit(Ω, 0.2, (1.0, -1.0))
+
+viz(Ωs.geometry)
+viz!(Ωt.geometry, color = :gray90)
+Mke.current_figure()
We can visualize the domain of the "train" (or source) set Ωs in blue, and the domain of the "test" (or target) set Ωt in gray. We reserved 20% of the samples to Ωs and 80% to Ωt. Internally, this geospatial geosplit function is implemented in terms of efficient geospatial partitions.
Let's define our geostatistical learning model to predict the crop type based on the four satellite bands. We will use the DecisionTreeClassifier model, which is suitable for the task we want to perform. Any model from the StatsLeanModels.jl model is supported, including all models from ScikitLearn.jl:
feats = [:band1, :band2, :band3, :band4]
+label = :crop
+
+model = DecisionTreeClassifier()DecisionTreeClassifier
+max_depth: -1
+min_samples_leaf: 1
+min_samples_split: 2
+min_purity_increase: 0.0
+pruning_purity_threshold: 1.0
+n_subfeatures: 0
+classes: nothing
+root: nothingWe will fit the model in Ωs where the features and labels are available and predict in Ωt where the features are available. The Learn transform automatically fits the model to the data:
learn = Learn(Ωs, model, feats => label)Learn transform
+├─ model: FittedModel{DecisionTreeClassifier}
+└─ input: [:band1, :band2, :band3, :band4]The transform can be called with new data to generate predictions:
Ω̂t = learn(Ωt)| crop | +geometry | +
|---|---|
| Categorical | +Point | +
| [NoUnits] | +🖈 Cartesian{NoDatum} | +
| 1 | +(x: 746.0 m, y: 717.0 m) | +
| 2 | +(x: 953.0 m, y: 291.0 m) | +
| 2 | +(x: 1050.0 m, y: 796.0 m) | +
| 2 | +(x: 787.0 m, y: 921.0 m) | +
| 2 | +(x: 947.0 m, y: 758.0 m) | +
| 2 | +(x: 949.0 m, y: 769.0 m) | +
| 2 | +(x: 773.0 m, y: 899.0 m) | +
| 2 | +(x: 851.0 m, y: 895.0 m) | +
| 3 | +(x: 1018.0 m, y: 748.0 m) | +
| 1 | +(x: 427.0 m, y: 743.0 m) | +
| ⋮ | +⋮ | +
We note that the prediction of a geostatistical learning model is a geospatial data object, and we can inspect it with the same methods already described above. This also means that we can visualize the prediction directly, side by side with the true label in this synthetic example:
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
With this example we conclude the basic workflow. To get familiar with other features of the project, please check the the reference guide.
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
Base.rand — Methodrand([rng], process::FieldProcess, domain, data, [nreals], [method]; [parameters])Generate one or nreals realizations of the field process with method over the domain with data and optional paramaters. Optionally, specify the random number generator rng and global parameters.
The data can be a geotable, a pair, or an iterable of pairs of the form var => T, where var is a symbol or string with the variable name and T is the corresponding data type.
Parameters
workers - Worker processes (default to workers())threads - Number of threads (default to cpucores())verbose - Show progress and other information (default to true)async - Return to master process immediately (default to false) Examples
julia> rand(process, domain, geotable, 3)
+julia> rand(process, domain, :z => Float64)
+julia> rand(process, domain, "z" => Float64)
+julia> rand(process, domain, [:a => Float64, :b => Float64])
+julia> rand(process, domain, ["a" => Float64, "b" => Float64])Realizations are stored in an Ensemble as illustrated in the following example:
GeoStatsProcesses.Ensemble — TypeEnsembleAn ensemble of geostatistical realizations from a geostatistical process.
# domain of interest
+grid = CartesianGrid(100, 100)
+
+# Gaussian process
+proc = GaussianProcess(GaussianVariogram(range=30.0))
+
+# unconditional simulation
+real = rand(proc, grid, [:Z => Float64], 100)2D Ensemble
+ domain: 100×100 CartesianGrid
+ variables: Z
+ N° reals: 100We can visualize the first two realizations:
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
the mean and variance:
m, v = mean(real), var(real)
+
+fig = Mke.Figure(size = (800, 400))
+viz(fig[1,1], m.geometry, color = m.Z)
+viz(fig[1,2], v.geometry, color = v.Z)
+fig
or the 25th and 75th percentiles:
q25 = quantile(real, 0.25)
+q75 = quantile(real, 0.75)
+
+fig = Mke.Figure(size = (800, 400))
+viz(fig[1,1], q25.geometry, color = q25.Z)
+viz(fig[1,2], q75.geometry, color = q75.Z)
+fig
The cummulative distribution function is obtained likewise:
p25 = cdf(real, 0.25)
+p75 = cdf(real, 0.75)
+
+fig = Mke.Figure(size = (800, 400))
+viz(fig[1,1], p25.geometry, color = p25.Z)
+viz(fig[1,2], p75.geometry, color = p75.Z)
+fig
All field processes can generate realizations in parallel using multiple Julia processes. Doing so requires using the Distributed standard library, like in the following example:
using Distributed
+
+# request additional processes
+addprocs(3)
+
+# load code on every single process
+@everywhere using GeoStats
+
+# setup simulation
+grid = CartesianGrid(100, 100)
+proc = GaussianProcess(GaussianVariogram(range=30.0))
+
+# perform simulation on all worker processes
+real = rand(proc, grid, [:Z => Float64], 3, workers = workers())Please consult The ultimate guide to distributed computing in Julia.
The following processes are shipped with the framework.
GeoStatsProcesses.GaussianProcess — TypeGaussianProcess(variogram=GaussianVariogram(), mean=0.0)Gaussian process with given theoretical variogram and global mean.
GeoStatsProcesses.LUMethod — TypeLUMethod(; [paramaters])The LU Gaussian process method introduced by Alabert 1987. The full covariance matrix is built to include all locations of the process domain, and samples from the multivariate Gaussian are drawn via LU factorization.
Parameters
factorization - Factorization method (default to cholesky)correlation - Correlation coefficient between two covariates (default to 0)init - Data initialization method (default to NearestInit())References
Alabert 1987. The practice of fast conditional simulations through the LU decomposition of the covariance matrix
Oliver 2003. Gaussian cosimulation: modeling of the cross-covariance
Notes
GeoStatsProcesses.FFTMethod — TypeFFTMethod(; [paramaters])The FFT Gaussian process method introduced by Gutjahr 1997. The covariance function is perturbed in the frequency domain after a fast Fourier transform. White noise is added to the phase of the spectrum, and a realization is produced by an inverse Fourier transform.
Parameters
minneighbors - Minimum number of neighbors (default to 1)maxneighbors - Maximum number of neighbors (default to 10)neighborhood - Search neighborhood (default to nothing)distance - Distance used to find nearest neighbors (default to Euclidean())References
Gutjahr 1997. General joint conditional simulations using a fast Fourier transform method
Gómez-Hernández, J. & Srivastava, R. 2021. One Step at a Time: The Origins of Sequential Simulation and Beyond
Notes
The method is limited to simulations on Cartesian grids, and care must be taken to make sure that the correlation length is small enough compared to the grid size.
As a general rule of thumb, avoid correlation lengths greater than 1/3 of the grid.
The method is extremely fast, and can be used to generate large 3D realizations.
GeoStatsProcesses.SEQMethod — TypeSEQMethod(; [paramaters])The sequential process method introduced by Gomez-Hernandez 1993. It traverses all locations of the geospatial domain according to a path, approximates the conditional Gaussian distribution within a neighborhood using simple Kriging, and assigns a value to the center of the neighborhood by sampling from this distribution.
Parameters
path - Process path (default to LinearPath())minneighbors - Minimum number of neighbors (default to 1)maxneighbors - Maximum number of neighbors (default to 36)neighborhood - Search neighborhood (default to :range)distance - Distance used to find nearest neighbors (default to Euclidean())init - Data initialization method (default to NearestInit())For each location in the process path, a maximum number of neighbors maxneighbors is used to fit the conditional Gaussian distribution. The neighbors are searched according to a neighborhood.
The neighborhood can be a MetricBall, the symbol :range or nothing. The symbol :range is converted to MetricBall(range(γ)) where γ is the variogram of the Gaussian process. If neighborhood is nothing, the nearest neighbors are used without additional constraints.
References
Notes
# domain of interest
+grid = CartesianGrid(100, 100)
+
+# Gaussian process
+proc = GaussianProcess(GaussianVariogram(range=30.0))
+
+# unconditional simulation
+real = rand(proc, grid, [:Z => Float64], 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
GeoStatsProcesses.LindgrenProcess — TypeLindgrenProcess(range=1.0, sill=1.0; init=NearestInit())Lindgren process with given range (correlation length) and sill (total variance) as described in Lindgren 2011. Optionally, specify the data initialization method init.
The process relies relies on a discretization of the Laplace-Beltrami operator on meshes and is adequate for highly curved domains (e.g. surfaces).
References
# domain of interest
+mesh = simplexify(Sphere((0, 0, 0), 1))
+
+# Lindgren process
+proc = LindgrenProcess()
+
+# unconditional simulation
+real = rand(proc, mesh, [:Z => Float64], 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
The following processes are available in external packages.
GeoStatsProcesses.QuiltingProcess — TypeQuiltingProcess(trainimg, tilesize; [paramaters])Image quilting process with training image trainimg and tile size tilesize as described in Hoffimann et al. 2017.
Parameters
overlap - Overlap size (default to (1/6, 1/6, ..., 1/6))path - Process path (:raster (default), :dilation, or :random)inactive - Vector of inactive voxels (i.e. CartesianIndex) in the gridsoft - A pair (data,dataTI) of geospatial data objects (default to nothing)tol - Initial relaxation tolerance in (0,1] (default to 0.1)init - Data initialization method (default to NearestInit())References
Hoffimann et al 2017. Stochastic simulation by image quilting of process-based geological models
Hoffimann et al 2015. Geostatistical modeling of evolving landscapes by means of image quilting
using ImageQuilting
+using GeoStatsImages
+
+# domain of interest
+grid = CartesianGrid(200, 200)
+
+# quilting process
+img = geostatsimage("Strebelle")
+proc = QuiltingProcess(img, (62, 62))
+
+# unconditional simulation
+real = rand(proc, grid, [:facies => Int], 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
# domain of interest
+grid = CartesianGrid(200, 200)
+
+# quilting process
+img = geostatsimage("StoneWall")
+proc = QuiltingProcess(img, (13, 13))
+
+# unconditional simulation
+real = rand(proc, grid, [:Z => Int], 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
# domain of interest
+grid = domain(img)
+
+# 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, data, 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
Voxels marked with the special symbol NaN are treated as inactive. The algorithm will skip tiles that only contain inactive voxels to save computation and will generate realizations that are consistent with the mask. This is particularly useful with complex 3D models that have large inactive portions.
# domain of interest
+grid = domain(img)
+
+# skip circle at the center
+nx, ny = size(grid)
+r = 100; circle = []
+for i in 1:nx, j in 1:ny
+ if (i-nx÷2)^2 + (j-ny÷2)^2 < r^2
+ push!(circle, CartesianIndex(i, j))
+ end
+end
+
+# quilting process
+proc = QuiltingProcess(img, (62, 62), inactive = circle)
+
+# unconditional simulation
+real = rand(proc, grid, [:facies => Float64], 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
It is possible to incorporate auxiliary variables to guide the selection of patterns from the training image.
using ImageFiltering
+
+# 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
+
+# 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), [:Z => Float64], 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
GeoStatsProcesses.TuringProcess — TypeTuringProcess(; [paramaters])Turing process as described in Turing 1952.
Parameters
params - basic parameters (default to nothing)blur - blur algorithm (default to nothing)edge - edge condition (default to nothing)iter - number of iterations (default to 100)References
using TuringPatterns
+
+# domain of interest
+grid = CartesianGrid(200, 200)
+
+# unconditional simulation
+real = rand(TuringProcess(), grid, [:z => Float64], 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
GeoStatsProcesses.StrataProcess — TypeStrataProcess(environment; [paramaters])Strata process with given geological environment as described in Hoffimann 2018.
Parameters
state - Initial geological statestack - Stacking scheme (:erosional or :depositional)nepochs - Number of epochs (default to 10)References
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, [:z => Float64], 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
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
Base.rand — Methodrand([rng], process::PointProcess, geometry, [nreals])
+rand([rng], process::PointProcess, domain, [nreals])Generate one or nreals realizations of the point process inside geometry or domain. Optionally specify the random number generator rng.
# geometry of interest
+sphere = Sphere((0, 0, 0), 1)
+
+# homogeneous Poisson process
+proc = PoissonProcess(5.0)
+
+# sample two point patterns
+pset = rand(proc, sphere, 2)
+
+fig = Mke.Figure(size = (800, 400))
+viz(fig[1,1], sphere)
+viz!(fig[1,1], pset[1], color = :black)
+viz(fig[1,2], sphere)
+viz!(fig[1,2], pset[2], color = :black)
+fig
GeoStatsProcesses.ishomogeneous — Functionishomogeneous(process::PointProcess)Tells whether or not the spatial point process process is homogeneous.
GeoStatsProcesses.BinomialProcess — TypeBinomialProcess(n)A Binomial point process with n points.
# 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
GeoStatsProcesses.PoissonProcess — TypePoissonProcess(λ)A Poisson process with intensity λ. For a homogeneous process, define λ as a constant real value, while for an inhomogeneous process, define λ as a function or vector of values. If λ is a vector, it is assumed that the process is associated with a Domain with the same number of elements as λ.
# geometry of interest
+box = Box((0, 0), (100, 100))
+
+# intensity function
+λ(p) = sum(to(p))^2 / 10000
+
+# homogeneous process
+proc₁ = PoissonProcess(0.5)
+
+# inhomogeneous process
+proc₂ = PoissonProcess(λ)
+
+# 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
GeoStatsProcesses.InhibitionProcess — TypeInhibitionProcess(δ)An inhibition point process with minimum distance δ.
# 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
GeoStatsProcesses.ClusterProcess — TypeClusterProcess(proc, ofun)A cluster process with parent process proc and offsprings generated with ofun. It is a function that takes a parent point and returns a point pattern from another point process.
ClusterProcess(proc, offs, gfun)Alternatively, specify the parent process proc, the offspring process offs and the geometry function gfun. It is a function that takes a parent point and returns a geometry or domain for the offspring process.
# geometry of interest
+box = Box((0, 0), (5, 5))
+
+# Matérn process
+proc₁ = ClusterProcess(
+ PoissonProcess(1),
+ PoissonProcess(1000),
+ p -> Ball(p, 0.2)
+)
+
+# inhomogeneous parent and offspring processes
+proc₂ = ClusterProcess(
+ PoissonProcess(p -> 0.1 * sum(to(p) .^ 2)),
+ p -> rand(PoissonProcess(x -> 5000 * sum((x - p).^2)), Ball(p, 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
Base.union — Methodp₁ ∪ p₂Return the union of point processes p₁ and p₂.
# 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
GeoStatsProcesses.thin — Functionthin(process, method)Thin spatial point process with thinning method.
GeoStatsProcesses.RandomThinning — TypeRandomThinning(p)Random thinning with retention probability p, which can be a constant probability value in [0,1] or a function mapping a point to a probability.
Examples
RandomThinning(0.5)
+RandomThinning(p -> sum(to(p)))# 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
# 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
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
The Julia ecosystem for geospatial modeling is maturing very quickly as the result of multiple initiatives such as JuliaEarth, JuliaClimate, and JuliaGeo. Each of these initiatives is associated with a different set of challenges that ultimatively determine the types of packages that are being developed in the corresponding GitHub organizations. In this section, we try to clarify what is available to first-time users of the language.
+
+
Originally created to host the GeoStats.jl framework, this initiative is primarily concerned with geospatial data science and geostatistical modeling. Due to the various applications in the subsurface of the Earth, most of our Julia packages were developed to work efficiently with both 2D and 3D geometries.
Unlike other initiatives, JuliaEarth is 100% Julia by design. This means that we do not rely on external libraries such as GDAL or Proj4 for geospatial work.
The most recent of the three initiatives, JuliaClimate has been created to address specific challenges in climate modeling. One of these challenges is access to climate data in Julia. Packages such as INMET.jl and CDSAPI.jl serve this purpose and are quite nice to work with.
Focused on bringing well-established external libraries to Julia, JuliaGeo provides packages that are widely used by geospatial communities from other programming languages. GDAL.jl, Proj4.jl and LibGEOS.jl are good examples of such packages.
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
We recommend the following educational resources.
Hoffimann, J. 2023. Geospatial Data Science with Julia - A fresh approach to data science with geospatial data and the Julia programming language.
Isaaks, E. and Srivastava, R. 1990. An Introduction to Applied Geostatistics - An introductory book on geostatistics that covers important concepts with very simple language.
Chilès, JP and Delfiner, P. 2012. Geostatistics: Modeling Spatial Uncertainty - An incredible book for those with good mathematical background.
Webster, R and Oliver, MA. 2007. Geostatistics for Environmental Scientists - A great book with good balance between theory and practice.
Journel, A. and Huijbregts, Ch. 2003. Mining Geostatistics - A great book with both theoretical and practical developments.
Júlio Hoffimann - Video lectures with the GeoStats.jl framework.
Edward Isaaks - Video lectures on variography, Kriging and related concepts.
Jef Caers - Video lectures on two-point and multiple-point methods.
UFMG 2024 [Portuguese] - Geociência de Dados na Mineração, UFMG 2024
JuliaEO 2024 [English] - Global Workshop on Earth Observation, AIRCentre 2024
UFMG 2023 [Portuguese] - Geociência de Dados na Mineração, UFMG 2023
JuliaEO 2023 [English] - Global Workshop on Earth Observation, AIRCentre 2023
CBMina 2021 [Portuguese] - Introução à Geoestatística, CBMina 2021
UFMG 2021 [Portuguese] - Introdução à Geoestatística, UFMG 2021
GaussianProcesses.jl - Gaussian process regression and Simple Kriging are essentially two names for the same concept. The derivation of Kriging estimators, however; does not require distributional assumptions. It is a beautiful coincidence that for multivariate Gaussian distributions, Simple Kriging gives the conditional expectation.
KernelFunctions.jl - Spatial structure can be represented in many different forms: covariance, variogram, correlogram, etc. Variograms are more general than covariance kernels according to the intrinsic stationary property. This means that there are variogram models with no covariance counterpart. Furthermore, empirical variograms can be easily estimated from the data (in various directions) with an efficient procedure. GeoStats.jl treats variograms as first-class objects.
Interpolations.jl - Kriging and spline interpolation have different purposes, yet these two methods are sometimes listed as competing alternatives. Kriging estimation is about minimizing variance (or estimation error), whereas spline interpolation is about deriving smooth estimators for computer visualization. Kriging is a generalization of splines in which one has the freedom to customize geospatial structure based on data. Besides the estimate itself, Kriging also provides the variance map as a function of point patterns.
ScikitLearn.jl - Traditional statistical learning relies on core assumptions that do not hold in geospatial settings (fixed support, i.i.d. samples, ...). Geostatistical learning has been introduced recently as an attempt to push the frontiers of statistical learning with geospatial data.
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
Below is a list of publications made possible with this project:
Kim et al. 2024. Coal Thickness Through Compositional Kriging: An Approach Based on Geostatistics
Doherty et al. 2024 A Method for Quantifying Uncertainty in Spatially Interpolated Meteorological Data with Application to Daily Maximum Air Temperature
Ko et al. 2024. Scaling Hawkes processes to one million COVID-19 cases
Hoffimann, J. 2023. Geospatial Data Science with Julia
Gruchalla et al. 2023. Reevaluating Contour Visualizations for Power Systems Data
Karsanina, M. & Gerke, K. 2022. Stochastic (re)constructions of non-stationary material structures: Using ensemble averaged correlation functions and non-uniform phase distributions
Sepúlveda et al. 2022. Evaluation of multivariate Gaussian transforms for geostatistical applications
Hoffimann et al. 2022. Modeling Geospatial Uncertainty of Geometallurgical Variables with Bayesian Models and Hilbert-Kriging
Hoffimann et al. 2021. Geostatistical Learning: Challenges and Opportunities
Asner et al. 2021. Abiotic and Human Drivers of Reef Habitat Complexity Throughout the Main Hawaiian Islands
Asner et al. 2020. Large-scale mapping of live corals to guide reef conservation
Hoffimann, J. & Zadrozny, B. 2019. Efficient Variography with Partition Variograms
Hoffimann et al. 2019. Morphodynamic Analysis and Statistical Synthesis of Geomorphic Data: Application to a Flume Experiment
Barfod et al. 2017. Hydrostratigraphic Modeling using Multiple-point Statistics and Airborne Transient Electromagnetic Methods
Hoffimann et al. 2017. Stochastic Simulation by Image Quilting of Process-based Geological Models
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
\n\n
","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"The Matheron's estimator of the empirical variogram is given by","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"widehatgamma_M(h) = frac12N(h) sum_(ij) in N(h) (z_i - z_j)^2","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"where N(h) = left(ij) mid x_i - x_j = hright is the set of pairs of locations at a distance h and N(h) is the cardinality of the set. Alternatively, the robust Cressie's estimator is given by","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"widehatgamma_C(h) = frac12fracleftfrac1N(h) sum_(ij) in N(h) z_i - z_j^12right^40457 + frac0494N(h) + frac0045N(h)^2","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"Both estimators are available and can be used with general distance functions in order to for example:","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"Model anisotropy (e.g. ellipsoid distance)\nPerform simulation on sphere (e.g. haversine distance)","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"Please see Distances.jl for a complete list of distance functions.","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"The high-performance estimation procedure implemented in the framework can consider all pairs of locations regardless of direction (ominidirectional) or a specified partition of the geospatial data (e.g. directional, planar).","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"Variograms can plotted with the following options:","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"varioplot","category":"page"},{"location":"variography/empirical/#GeoStatsFunctions.varioplot","page":"Empirical variograms","title":"GeoStatsFunctions.varioplot","text":"varioplot(γ; [options])\n\nPlot the variogram γ with given options.\n\nCommon variogram options:\n\ncolor - color of variogram\nsize - size of variogram\nmaxlag - maximum lag of variogram\n\nEmpirical variogram options:\n\npointsize - size of points of variogram\nshowtext - show text counts\ntextsize - size of text counts\nshowhist - show histogram\nhistcolor - color of histogram\n\nNotes\n\nThis function will only work in the presence of a Makie.jl backend via package extensions in Julia v1.9 or later versions of the language.\n\n\n\n\n\n","category":"function"},{"location":"variography/empirical/#(Omini)directional","page":"Empirical variograms","title":"(Omini)directional","text":"","category":"section"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"EmpiricalVariogram\nDirectionalVariogram\nPlanarVariogram","category":"page"},{"location":"variography/empirical/#GeoStatsFunctions.EmpiricalVariogram","page":"Empirical variograms","title":"GeoStatsFunctions.EmpiricalVariogram","text":"EmpiricalVariogram(data, var₁, var₂=var₁; [parameters])\n\nComputes the empirical (a.k.a. experimental) omnidirectional (cross-)variogram for variables var₁ and var₂ stored in geospatial data.\n\nParameters\n\nnlags - number of lags (default to 20)\nmaxlag - maximum lag in length units (default to 1/2 of minimum side of bounding box)\ndistance - custom distance function (default to Euclidean distance)\nestimator - variogram estimator (default to :matheron estimator)\nalgorithm - accumulation algorithm (default to :ball)\n\nAvailable estimators:\n\n:matheron - simple estimator based on squared differences\n:cressie - robust estimator based on 4th power of differences\n\nAvailable algorithms:\n\n:full - loop over all pairs of points in the data\n:ball - loop over all points inside maximum lag ball\n\nAll implemented algorithms produce the exact same result. The :ball algorithm is considerably faster when the maximum lag is much smaller than the bounding box of the domain of the data.\n\nSee also: DirectionalVariogram, PlanarVariogram, EmpiricalVarioplane.\n\nReferences\n\nChilès, JP and Delfiner, P. 2012. Geostatistics: Modeling Spatial Uncertainty\nWebster, R and Oliver, MA. 2007. Geostatistics for Environmental Scientists\nHoffimann, J and Zadrozny, B. 2019. Efficient variography with partition variograms\n\n\n\n\n\n","category":"type"},{"location":"variography/empirical/#GeoStatsFunctions.DirectionalVariogram","page":"Empirical variograms","title":"GeoStatsFunctions.DirectionalVariogram","text":"DirectionalVariogram(direction, data, var₁, var₂=var₁; dtol=1e-6u\"m\", [parameters])\n\nComputes the empirical (cross-)variogram for the variables var₁ and var₂ stored in geospatial data along a given direction with band tolerance dtol in length units.\n\nOptionally, forward parameters for the underlying EmpiricalVariogram.\n\n\n\n\n\n","category":"function"},{"location":"variography/empirical/#GeoStatsFunctions.PlanarVariogram","page":"Empirical variograms","title":"GeoStatsFunctions.PlanarVariogram","text":"PlanarVariogram(normal, data, var₁, var₂=var₁; ntol=1e-6u\"m\", [parameters])\n\nComputes the empirical (cross-)variogram for the variables var₁ and var₂ stored in geospatial data along a plane perpendicular to a normal direction with plane tolerance ntol in length units.\n\nOptionally, forward parameters for the underlying EmpiricalVariogram.\n\n\n\n\n\n","category":"function"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"Consider the following example image:","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"using GeoStatsImages\n\nimg = geostatsimage(\"Gaussian30x10\")\n\nimg |> viewer","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"We can estimate ominidirectional variograms, which consider pairs of points along all directions:","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"γ = EmpiricalVariogram(img, :Z, maxlag = 50.)\n\nvarioplot(γ)","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"directional variograms along a specific direction:","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"γₕ = DirectionalVariogram((1.,0.), img, :Z, maxlag = 50.)\nγᵥ = DirectionalVariogram((0.,1.), img, :Z, maxlag = 50.)\n\nvarioplot(γₕ, color = :maroon, histcolor = :maroon)\nvarioplot!(γᵥ)\nMke.current_figure()","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"or planar variograms over a specific plane:","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"γᵥ = PlanarVariogram((1.,0.), img, :Z, maxlag = 50.)\nγₕ = PlanarVariogram((0.,1.), img, :Z, maxlag = 50.)\n\nvarioplot(γₕ, color = :maroon, histcolor = :maroon)\nvarioplot!(γᵥ)\nMke.current_figure()","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"note: Note\nThe directional and planar variograms coincide in this example because planes are equal to lines in 2-dimensional space. These concepts are most useful in 3-dimensional space where we may be interested in comparing the horizontal planar range to the vertical directional range.","category":"page"},{"location":"variography/empirical/#Varioplanes","page":"Empirical variograms","title":"Varioplanes","text":"","category":"section"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"Variograms estimated along all directions in a given plane of reference are called varioplanes.","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"EmpiricalVarioplane","category":"page"},{"location":"variography/empirical/#GeoStatsFunctions.EmpiricalVarioplane","page":"Empirical variograms","title":"GeoStatsFunctions.EmpiricalVarioplane","text":"EmpiricalVarioplane(data, var₁, var₂=var₁;\n normal=Vec(0,0,1), nangs=50,\n ptol=0.5u\"m\", dtol=0.5u\"m\",\n [parameters])\n\nGiven a normal direction, estimate the (cross-)variogram of variables var₁ and var₂ along all directions in the corresponding plane of variation.\n\nOptionally, specify the tolerance ptol in length units for the plane partition, the tolerance dtol in length units for the direction partition, the number of angles nangs in the plane, and forward the parameters to the underlying EmpiricalVariogram.\n\n\n\n\n\n","category":"type"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"The varioplane is plotted on a polar axis for all lags and angles:","category":"page"},{"location":"variography/empirical/","page":"Empirical variograms","title":"Empirical variograms","text":"γ = EmpiricalVarioplane(img, :Z, maxlag = 50.)\n\nplaneplot(γ)","category":"page"},{"location":"variography/theoretical/#Theoretical-variograms","page":"Theoretical variograms","title":"Theoretical variograms","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"using GeoStats # hide\nimport CairoMakie as Mke # hide","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"We provide various theoretical variogram models from the literature, which can can be combined with ellipsoid distances to model geometric anisotropy and nested with scalars or matrix coefficients to express multivariate relations.","category":"page"},{"location":"variography/theoretical/#Models","page":"Theoretical variograms","title":"Models","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"In an intrinsic isotropic model, the variogram is only a function of the distance between any two points x_1x_2 in R^m:","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(x_1x_2) = gamma(x_1 - x_2) = gamma(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"Under the additional assumption of 2nd-order stationarity, the well-known covariance is directly related via gamma(h) = cov(0) - cov(h). This package implements a few commonly used as well as other more exotic variogram models. Most of these models share a set of default parameters (e.g. sill=1.0, range=1.0), which can be set with keyword arguments.","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"Functions are provided to query properties of variogram models programmatically:","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"isstationary(::Variogram)\nisisotropic(::Variogram)","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.isstationary-Tuple{Variogram}","page":"Theoretical variograms","title":"GeoStatsFunctions.isstationary","text":"isstationary(f)\n\nCheck if geostatistical function f possesses the 2nd-order stationary property.\n\n\n\n\n\n","category":"method"},{"location":"variography/theoretical/#Meshes.isisotropic-Tuple{Variogram}","page":"Theoretical variograms","title":"Meshes.isisotropic","text":"isisotropic(f)\n\nTells whether or not the geostatistical function f is isotropic.\n\n\n\n\n\n","category":"method"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"Corresponding covariance models are available for convenience. For example, the GaussianVariogram and the GaussianCovariance are two alternative representations of the same geostatistical behavior.","category":"page"},{"location":"variography/theoretical/#Gaussian","page":"Theoretical variograms","title":"Gaussian","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(h) = (s - n) left1 - expleft(-3left(frachrright)^2right)right + n cdot 1_(0infty)(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"GaussianVariogram","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.GaussianVariogram","page":"Theoretical variograms","title":"GeoStatsFunctions.GaussianVariogram","text":"GaussianVariogram(; range, sill, nugget)\nGaussianVariogram(ball; sill, nugget)\n\nA Gaussian variogram with range, sill and nugget. Optionally, use a custom metric ball.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"varioplot(GaussianVariogram())","category":"page"},{"location":"variography/theoretical/#Spherical","page":"Theoretical variograms","title":"Spherical","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(h) = (s - n) leftleft(frac32left(frachrright) + frac12left(frachrright)^3right) cdot 1_(0r)(h) + 1_rinfty)(h)right + n cdot 1_(0infty)(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"SphericalVariogram","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.SphericalVariogram","page":"Theoretical variograms","title":"GeoStatsFunctions.SphericalVariogram","text":"SphericalVariogram(; range, sill, nugget)\nSphericalVariogram(ball; sill, nugget)\n\nA spherical variogram with range, sill and nugget. Optionally, use a custom metric ball.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"varioplot(SphericalVariogram())","category":"page"},{"location":"variography/theoretical/#Exponential","page":"Theoretical variograms","title":"Exponential","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(h) = (s - n) left1 - expleft(-3left(frachrright)right)right + n cdot 1_(0infty)(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"ExponentialVariogram","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.ExponentialVariogram","page":"Theoretical variograms","title":"GeoStatsFunctions.ExponentialVariogram","text":"ExponentialVariogram(; range, sill, nugget)\nExponentialVariogram(ball; sill, nugget)\n\nA exponential variogram with range, sill and nugget. Optionally, use a custom metric ball.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"varioplot(ExponentialVariogram())","category":"page"},{"location":"variography/theoretical/#Matern","page":"Theoretical variograms","title":"Matern","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(h) = (s - n) left1 - frac2^1-nuGamma(nu) left(sqrt2nu 3left(frachrright)right)^nu K_nuleft(sqrt2nu 3left(frachrright)right)right + n cdot 1_(0infty)(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"MaternVariogram","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.MaternVariogram","page":"Theoretical variograms","title":"GeoStatsFunctions.MaternVariogram","text":"MaternVariogram(; range, sill, nugget, order)\nMaternVariogram(ball; sill, nugget, order)\n\nA Matérn variogram with range, sill and nugget. Optionally, use a custom metric ball and order of the Bessel function.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"varioplot(MaternVariogram())","category":"page"},{"location":"variography/theoretical/#Cubic","page":"Theoretical variograms","title":"Cubic","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(h) = (s - n) leftleft(7left(frachrright)^2 - frac354left(frachrright)^3 + frac72left(frachrright)^5 - frac34left(frachrright)^7right) cdot 1_(0r)(h) + 1_rinfty)(h)right + n cdot 1_(0infty)(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"CubicVariogram","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.CubicVariogram","page":"Theoretical variograms","title":"GeoStatsFunctions.CubicVariogram","text":"CubicVariogram(; range, sill, nugget)\nCubicVariogram(ball; sill, nugget)\n\nA cubic variogram with range, sill and nugget. Optionally, use a custom metric ball.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"varioplot(CubicVariogram())","category":"page"},{"location":"variography/theoretical/#Pentaspherical","page":"Theoretical variograms","title":"Pentaspherical","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(h) = (s - n) leftleft(frac158left(frachrright) - frac54left(frachrright)^3 + frac38left(frachrright)^5right) cdot 1_(0r)(h) + 1_rinfty)(h)right + n cdot 1_(0infty)(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"PentasphericalVariogram","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.PentasphericalVariogram","page":"Theoretical variograms","title":"GeoStatsFunctions.PentasphericalVariogram","text":"PentasphericalVariogram(; range, sill, nugget)\nPentasphericalVariogram(ball; sill, nugget)\n\nA pentaspherical variogram with range, sill and nugget. Optionally, use a custom metric ball.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"varioplot(PentasphericalVariogram())","category":"page"},{"location":"variography/theoretical/#Sine-hole","page":"Theoretical variograms","title":"Sine hole","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(h) = (s - n) left1 - fracsin(pi h r)pi h rright + n cdot 1_(0infty)(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"SineHoleVariogram","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.SineHoleVariogram","page":"Theoretical variograms","title":"GeoStatsFunctions.SineHoleVariogram","text":"SineHoleVariogram(; range, sill, nugget)\nSineHoleVariogram(ball; sill, nugget)\n\nA sine hole variogram with range, sill and nugget. Optionally, use a custom metric ball.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"varioplot(SineHoleVariogram())","category":"page"},{"location":"variography/theoretical/#Circular","page":"Theoretical variograms","title":"Circular","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(h) = (s - n) leftleft(1 - frac2pi cos^-1left(frachrright) + frac2hpi r sqrt1 - frach^2r^2 right) cdot 1_(0r)(h) + 1_rinfty)(h)right + n cdot 1_(0infty)(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"CircularVariogram","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.CircularVariogram","page":"Theoretical variograms","title":"GeoStatsFunctions.CircularVariogram","text":"CircularVariogram(; range, sill, nugget)\nCircularVariogram(ball; sill, nugget)\n\nA circular variogram with range, sill and nugget. Optionally, use a custom metric ball.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"varioplot(CircularVariogram())","category":"page"},{"location":"variography/theoretical/#Power","page":"Theoretical variograms","title":"Power","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(h) = sh^a + n cdot 1_(0infty)(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"PowerVariogram","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.PowerVariogram","page":"Theoretical variograms","title":"GeoStatsFunctions.PowerVariogram","text":"PowerVariogram(; scaling, exponent, nugget)\n\nA power variogram with scaling, exponent and nugget.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"varioplot(PowerVariogram())","category":"page"},{"location":"variography/theoretical/#Nugget","page":"Theoretical variograms","title":"Nugget","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"gamma(h) = n cdot 1_(0infty)(h)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"NuggetEffect","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.NuggetEffect","page":"Theoretical variograms","title":"GeoStatsFunctions.NuggetEffect","text":"NuggetEffect(; nugget)\n\nA pure nugget effect variogram with nugget.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"varioplot(NuggetEffect(1.0))","category":"page"},{"location":"variography/theoretical/#Anisotropy","page":"Theoretical variograms","title":"Anisotropy","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"Anisotropic models are easily obtained by defining an ellipsoid metric in place of the default Euclidean metric as shown in the following example. First, we create an ellipsoid that specifies the ranges and angles of rotation:","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"ellipsoid = MetricBall((3.0, 2.0, 1.0), RotZXZ(0.0, 0.0, 0.0))","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"All rotations from Rotations.jl are supported as well as the following additional rotations from commercial geostatistical software:","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"MinesightAngles\nDatamineAngles\nVulcanAngles\nGslibAngles","category":"page"},{"location":"variography/theoretical/#GeoStatsBase.MinesightAngles","page":"Theoretical variograms","title":"GeoStatsBase.MinesightAngles","text":"MinesightAngles(θ₁, θ₂, θ₃)\n\nMineSight z-x'-y'' intrinsic rotation convention following the right-hand rule. All angles are in degrees and the sign convention is CW, CCW, CW positive.\n\nThe first rotation θ₁ is a horizontal rotation around the Z-axis, with positive being clockwise. The second rotation θ₂ is a rotation around the new X-axis, which moves the Y-axis into the desired position, with positive direction of rotation is up. The third rotation θ₃ is a rotation around the new Y-axis, which moves the X-axis into the desired position, with positive direction of rotation is up.\n\nReferences\n\nSanchez, J. MINESIGHT® TUTORIALS\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/#GeoStatsBase.DatamineAngles","page":"Theoretical variograms","title":"GeoStatsBase.DatamineAngles","text":"DatamineAngles(θ₁, θ₂, θ₃)\n\nDatamine ZXY rotation convention following the left-hand rule. All angles are in degrees and the signal convention is CW, CW, CW positive. Y is the principal axis.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/#GeoStatsBase.VulcanAngles","page":"Theoretical variograms","title":"GeoStatsBase.VulcanAngles","text":"VulcanAngles(θ₁, θ₂, θ₃)\n\nVulcan ZYX rotation convention following the right-hand rule. All angles are in degrees and the signal convention is CW, CCW, CW positive. X is the principal axis.\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/#GeoStatsBase.GslibAngles","page":"Theoretical variograms","title":"GeoStatsBase.GslibAngles","text":"GslibAngles(θ₁, θ₂, θ₃)\n\nGSLIB z-x'-y'' intrinsic rotation convention following the left-hand rule. All angles are in degrees and the sign convention is CW, CCW, CW positive. Y is the principal axis.\n\nThe first rotation θ₁ is a rotation around the Z-axis, this is also called the azimuth. The second rotation θ₂ is a rotation around the new X-axis, this is also called the dip. The third rotation θ₃ is a rotation around the new Y-axis, this is also called the tilt. \n\nReferences\n\nDeutsch, 2015. The Angle Specification for GSLIB Software\n\n\n\n\n\n","category":"type"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"We pass the ellipsoid as the first argument to the variogram model instead of specifying a single range with a keyword argument:","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"GaussianVariogram(ellipsoid, sill = 2.0)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"To illustrate the concept, consider the following 2D data set:","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"using Random # hide\nRandom.seed!(2000) # hide\n\npoints = rand(Point, 50, crs=Cartesian2D) |> Scale(100)\ngeotable = georef((Z=rand(50),), points)\n\nviz(geotable.geometry, color = geotable.Z)","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"We interpolate the data with different ellipsoids by varying the angle of rotation from 0 to 2pi clockwise:","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"θs = range(0.0, step = π/4, stop = 2π - π/4)\n\n# initialize figure\nfig = Mke.Figure(size = (800, 1600))\n\n# helper function to position subfigures\npos = i -> CartesianIndices((4, 2))[i].I\n\n# Kriging with different angles\nfor (i, θ) in enumerate(θs)\n # ellipsoid rotated clockwise by angle θ\n e = MetricBall((20.,5.), Angle2d(θ))\n\n # anisotropic variogram model\n γ = GaussianVariogram(e)\n\n # domain of interpolation\n grid = CartesianGrid(100, 100)\n\n # perform interpolation\n sol = geotable |> Interpolate(grid, Kriging(γ))\n\n # visualize solution at subplot i\n viz(fig[pos(i)...], sol.geometry, color = sol.Z,\n axis = (title = \"θ = $(rad2deg(θ))ᵒ\",)\n )\nend\n\nfig","category":"page"},{"location":"variography/theoretical/#Nesting","page":"Theoretical variograms","title":"Nesting","text":"","category":"section"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"A nested variogram model gamma(h) = c_1cdotgamma_1(h) + c_2cdotgamma_2(h) + cdots + c_ncdotgamma_n(h) can be constructed from multiple variogram models, including matrix coefficients. The individual structures can be recovered in canonical form with the structures function:","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"structures","category":"page"},{"location":"variography/theoretical/#GeoStatsFunctions.structures","page":"Theoretical variograms","title":"GeoStatsFunctions.structures","text":"structures(γ)\n\nReturn the individual structures of a (possibly nested) variogram as a tuple. The structures are the total nugget cₒ, and the coefficients (or contributions) c[i] for the remaining non-trivial structures g[i] after normalization (i.e. sill=1, nugget=0).\n\nExamples\n\nγ₁ = GaussianVariogram(nugget=1, sill=2)\nγ₂ = SphericalVariogram(nugget=2, sill=3)\n\nγ = 2γ₁ + 3γ₂\n\ncₒ, c, g = structures(γ)\n\n\n\n\n\n","category":"function"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"γ₁ = GaussianVariogram(nugget=1, sill=2)\nγ₂ = ExponentialVariogram(nugget=2, sill=3)\n\n# nested model with matrix coefficients\nγ = [1.0 0.0; 0.0 1.0] * γ₁ + [2.0 0.5; 0.5 3.0] * γ₂\n\n# structures in canonical form\ncₒ, c, g = structures(γ)\n\ncₒ # matrix nugget","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"c # matrix coefficients","category":"page"},{"location":"variography/theoretical/","page":"Theoretical variograms","title":"Theoretical variograms","text":"g # normalized structures","category":"page"},{"location":"contributing/#Contributing","page":"Contributing","title":"Contributing","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"First off, thank you for considering contributing to GeoStats.jl. It’s people like you that make this project so much fun. Below are a few suggestions to speed up the collaboration process:","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Please be polite, we are here to help and learn from each other.\nTry to explain your contribution with simple language.\nReferences to textbooks and papers are always welcome.\nFollow the coding standards in the source.","category":"page"},{"location":"contributing/#How-to-start-contributing?","page":"Contributing","title":"How to start contributing?","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Contributing to an open-source project for the very first time can be a very daunting task. To make the process easier and more GitHub-beginner-friendly, the community has written an article about how to start contributing to open-source and overcome the mental and technical barriers that come associated with it. The article will also take you through the steps required to make your first contribution in detail.","category":"page"},{"location":"contributing/#Reporting-issues","page":"Contributing","title":"Reporting issues","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"If you are experiencing issues or have discovered a bug, please report it on GitHub. To make the resolution process easier, please include the version of Julia and GeoStats.jl in your writeup. These can be found with two commands:","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"julia> versioninfo()\njulia> using Pkg; Pkg.status()","category":"page"},{"location":"contributing/#Feature-requests","page":"Contributing","title":"Feature requests","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"If you have suggestions of improvement or algorithms that you would like to see implemented, please open an issue on GitHub. Suggestions as well as feature requests are very welcome.","category":"page"},{"location":"contributing/#Code-contribution","page":"Contributing","title":"Code contribution","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"If you have code that you would like to contribute that is awesome! Please open an issue or reach out in our community channel before you create the pull request on GitHub so that we make sure your idea is aligned with the project goals.","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"After your idea is revised by project maintainers, you implement it online on Github or offline on your machine.","category":"page"},{"location":"contributing/#Online-changes","page":"Contributing","title":"Online changes","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"If the changes to the code are minimal, we recommend pressing . on the keyboard on any file in the GitHub repository of interest. This will open the VSCode editor on the web browser where you can implement the changes, commit them and submit a pull request.","category":"page"},{"location":"contributing/#Offline-changes","page":"Contributing","title":"Offline changes","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"If the changes require additional investigation and tests, please get the development version of the project by typing the following in the package manager:","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"] activate @geo","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"This will create a fresh environment called @geo where you can edit the project modules without effects on your global user environment. Next, go ahead and ask the package manager to develop the package of interest (e.g. GeoStatsFunctions.jl):","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"] dev GeoStatsFunctions","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"You can modify the source code that was cloned in the .julia/dev folder and submit a pull request on GitHub later.","category":"page"},{"location":"links/#Index","page":"Index","title":"Index","text":"","category":"section"},{"location":"links/","page":"Index","title":"Index","text":"Below is the list of types and functions mentioned in the documentation.","category":"page"},{"location":"links/#Types","page":"Index","title":"Types","text":"","category":"section"},{"location":"links/","page":"Index","title":"Index","text":"Order = [:type]","category":"page"},{"location":"links/#Functions","page":"Index","title":"Functions","text":"","category":"section"},{"location":"links/","page":"Index","title":"Index","text":"Order = [:function]","category":"page"},{"location":"data/#Geospatial-data","page":"Data","title":"Geospatial data","text":"","category":"section"},{"location":"data/#Overview","page":"Data","title":"Overview","text":"","category":"section"},{"location":"data/","page":"Data","title":"Data","text":"Given a table or array containing data, we can georeference these objects onto a geospatial domain with the georef function. The result of the georef function is a GeoTable. If you would like to learn this concept in more depth, check out the recording of our JuliaEO 2024 workshop:","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"\n\n
","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"georef","category":"page"},{"location":"data/#GeoTables.georef","page":"Data","title":"GeoTables.georef","text":"georef(table, domain)\n\nGeoreference table on domain from Meshes.jl\n\nExamples\n\njulia> georef((a=rand(100), b=rand(100)), CartesianGrid(10, 10))\n\n\n\n\n\ngeoref(table, geoms)\n\nGeoreference table on vector of geometries geoms from Meshes.jl\n\nExamples\n\njulia> georef((a=rand(10), b=rand(10)), rand(Point, 10))\n\n\n\n\n\ngeoref(table, coords; [crs], [lenunit])\n\nGeoreference table using coordinates coords of points.\n\nOptionally, specify the coordinate reference system crs, which is set by default based on heuristics, and the lenunit (default to meters for unitless values) that can only be used in CRS types that allow this flexibility. Any CRS or EPSG/ESRI code from CoordRefSystems.jl is supported.\n\nExamples\n\njulia> georef((a=[1, 2, 3], b=[4, 5, 6], [(0, 0), (1, 1), (2, 2)])\njulia> georef((a=[1, 2, 3], b=[4, 5, 6], [(0, 0), (1, 1), (2, 2)], crs=LatLon)\njulia> georef((a=[1, 2, 3], b=[4, 5, 6], [(0, 0), (1, 1), (2, 2)], crs=EPSG{4326})\njulia> georef((a=[1, 2, 3], b=[4, 5, 6], [(0, 0), (1, 1), (2, 2)], lenunit=u\"cm\")\n\n\n\n\n\ngeoref(table, names; [crs], [lenunit])\n\nGeoreference table using coordinates of points stored in column names.\n\nOptionally, specify the coordinate reference system crs, which is set by default based on heuristics, and the lenunit (default to meters for unitless values) that can only be used in CRS types that allow this flexibility. Any CRS or EPSG/ESRI code from CoordRefSystems.jl is supported.\n\nExamples\n\ngeoref((a=rand(10), x=rand(10), y=rand(10)), (\"x\", \"y\"))\ngeoref((a=rand(10), x=rand(10), y=rand(10)), (\"x\", \"y\"), crs=LatLon)\ngeoref((a=rand(10), x=rand(10), y=rand(10)), (\"x\", \"y\"), crs=EPSG{4326})\ngeoref((a=rand(10), x=rand(10), y=rand(10)), (\"x\", \"y\"), lenunit=u\"cm\")\n\n\n\n\n\ngeoref(tuple)\n\nGeoreference a named tuple on CartesianGrid(dims), with dims obtained from the arrays stored in the tuple.\n\nExamples\n\njulia> georef((a=rand(10, 10), b=rand(10, 10))) # 2D grid\njulia> georef((a=rand(10, 10, 10), b=rand(10, 10, 10))) # 3D grid\n\n\n\n\n\n","category":"function"},{"location":"data/","page":"Data","title":"Data","text":"The functions values and domain can be used to retrieve the table of attributes and the underlying geospatial domain:","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"values\ndomain","category":"page"},{"location":"data/#Base.values","page":"Data","title":"Base.values","text":"values(geotable, [rank])\n\nReturn the values of geotable for a given rank as a table.\n\nThe rank is a non-negative integer that specifies the parametric dimension of the geometries of interest:\n\n0 - points\n1 - segments\n2 - triangles, quadrangles, ...\n3 - tetrahedrons, hexahedrons, ...\n\nIf the rank is not specified, it is assumed to be the rank of the elements of the domain.\n\n\n\n\n\n","category":"function"},{"location":"data/#GeoTables.domain","page":"Data","title":"GeoTables.domain","text":"domain(geotable)\n\nReturn underlying domain of the geotable.\n\n\n\n\n\n","category":"function"},{"location":"data/#Examples","page":"Data","title":"Examples","text":"","category":"section"},{"location":"data/","page":"Data","title":"Data","text":"using GeoStats\nimport CairoMakie as Mke\n\n# helper function for plotting two\n# variables named T and P side by side\nfunction plot(data)\n fig = Mke.Figure(size = (800, 400))\n viz(fig[1,1], data.geometry, color = data.T)\n viz(fig[1,2], data.geometry, color = data.P)\n fig\nend","category":"page"},{"location":"data/#Tables","page":"Data","title":"Tables","text":"","category":"section"},{"location":"data/","page":"Data","title":"Data","text":"Consider a table (e.g. DataFrame) with 25 samples of temperature T and pressure P:","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"using DataFrames\n\ntable = DataFrame(T=rand(25), P=rand(25))","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"We can georeference this table based on a given set of points:","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"georef(table, rand(Point, 25)) |> plot","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"or alternatively, georeference it on a 5x5 regular grid (5x5 = 25 samples):","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"georef(table, CartesianGrid(5, 5)) |> plot","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"Another common pattern in geospatial data is when the coordinates of the samples are already part of the table as columns. In this case, we can specify the column names:","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"table = DataFrame(T=rand(25), P=rand(25), X=rand(25), Y=rand(25), Z=rand(25))\n\ngeoref(table, (\"X\", \"Y\", \"Z\")) |> plot","category":"page"},{"location":"data/#Arrays","page":"Data","title":"Arrays","text":"","category":"section"},{"location":"data/","page":"Data","title":"Data","text":"Consider arrays (e.g. images) with data for various geospatial variables. We can georeference these arrays using a named tuple, and the framework will understand that the shape of the arrays should be preserved in a CartesianGrid:","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"T, P = rand(5, 5), rand(5, 5)\n\ngeoref((T=T, P=P)) |> plot","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"Alternatively, we can interpret the entries of the named tuple as columns in a table:","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"georef((T=vec(T), P=vec(P)), rand(Point, 25)) |> plot","category":"page"},{"location":"data/#Files","page":"Data","title":"Files","text":"","category":"section"},{"location":"data/","page":"Data","title":"Data","text":"The GeoIO.jl package can be used to load/save geospatial data from/to various file formats:","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"GeoIO.load\nGeoIO.save\nGeoIO.formats","category":"page"},{"location":"data/#GeoIO.load","page":"Data","title":"GeoIO.load","text":"load(fname, repair=true, layer=0, lenunit=m, kwargs...)\n\nLoad geospatial table from file fname stored in any format.\n\nVarious repairs are performed on the stored geometries by default, including fixes of orientation in rings of polygons, removal of zero-area triangles, etc.\n\nSome of the repairs can be expensive on large data sets. In that case, we recommend setting repair=false. Custom repairs can be performed with the Repair transform from Meshes.jl.\n\nOptionally, specify the layer to read within the file, and the length unit lenunit of the coordinates when the format does not include units in its specification. Other kwargs are forwarded to the backend packages.\n\nPlease use the formats function to list all supported file formats.\n\nOptions\n\nOFF\n\ndefaultcolor: default color of the geometries if the file does not have this data (default to RGBA(0.666, 0.666, 0.666, 0.666));\n\nCSV\n\ncoords: names of the columns with point coordinates (required option);\nOther options are passed to CSV.File, see the CSV.jl documentation for more details;\n\nVTK formats (.vtu, .vtp, .vtr, .vts, .vti)\n\nmask: name of the boolean column that encodes the indices of a grid view (default to :MASK). If the column does not exist in the file, the full grid is returned;\n\nCommon Data Model formats (NetCDF, GRIB)\n\nx: name of the column with x coordinates (default to \"x\", \"X\", \"lon\", or \"longitude\");\ny: name of the column with y coordinates (default to \"y\", \"Y\", \"lat\", or \"latitude\");\nz: name of the column with z coordinates (default to \"z\", \"Z\", \"depth\", or \"height\");\nt: name of the column with time measurements (default to \"t\", \"time\", or \"TIME\");\n\nGeoJSON\n\nnumbertype: number type of geometry coordinates (default to Float64)\nOther options are passed to GeoJSON.read, see the GeoJSON.jl documentation for more details;\n\nGSLIB\n\nOther options are passed to GslibIO.load, see the GslibIO.jl documentation for more details;\n\nShapefile\n\nOther options are passed to Shapefile.read, see the Shapefile.jl documentation for more details;\n\nGeoParquet\n\nOther options are passed to GeoParquet.read, see the GeoParquet.jl documentation for more details;\n\nGeoTIFF, GeoPackage, KML\n\nOther options are passed to ArchGDAL.read, see the ArchGDAL.jl documentation for more details;\n\nExamples\n\n# load coordinates of geojson file as Float64 (default)\nGeoIO.load(\"file.geojson\")\n# load coordinates of geojson file as Float32\nGeoIO.load(\"file.geojson\", numbertype=Float32)\n\n\n\n\n\n","category":"function"},{"location":"data/#GeoIO.save","page":"Data","title":"GeoIO.save","text":"save(fname, geotable; kwargs...)\n\nSave geotable to file fname of given format based on the file extension.\n\nOther kwargs are forwarded to the backend packages.\n\nPlease use the formats function to list all supported file formats.\n\nOptions\n\nOFF\n\ncolor: name of the column with geometry colors, if nothing the geometries will be saved without colors (default to nothing);\n\nMSH\n\nvcolumn: name of the column in vertex table with node data, if nothing the geometries will be saved without node data (default to nothing);\necolumn: name of the column in element table with element data, if nothing the geometries will be saved without element data (default to nothing);\n\nSTL\n\nascii: defines whether the file will be saved in ASCII format, otherwise Binary format will be used (default to false);\n\nCSV\n\ncoords: names of the columns where the point coordinates will be saved (default to \"x\", \"y\", \"z\");\nfloatformat: C-style format string for float values (default to no formatting);\nOther options are passed to CSV.write, see the CSV.jl documentation for more details;\n\nNetCDF\n\nx: name of the column where the coordinate x will be saved (default to CRS coordinate name);\ny: name of the column where the coordinate y will be saved (default to CRS coordinate name);\nz: name of the column where the coordinate z will be saved (default to CRS coordinate name);\nt: name of the column where the time measurements will be saved (default to \"t\");\n\nGeoTIFF\n\noptions: list with options that will be passed to GDAL;\n\nGeoPackage\n\nlayername: name of the layer where the data will be saved (default to \"data\");\noptions: dictionary with options that will be passed to GDAL;\n\nGSLIB\n\nOther options are passed to GslibIO.save, see the GslibIO.jl documentation for more details;\n\nShapefile\n\nOther options are passed to Shapefile.write, see the Shapefile.jl documentation for more details;\n\nGeoJSON\n\nOther options are passed to GeoJSON.write, see the GeoJSON.jl documentation for more details;\n\nGeoParquet\n\nOther options are passed to GeoParquet.write, see the GeoParquet.jl documentation for more details;\n\nExamples\n\n# overwrite an existing shapefile\nGeoIO.save(\"file.shp\", force = true)\n\n\n\n\n\n","category":"function"},{"location":"data/#GeoIO.formats","page":"Data","title":"GeoIO.formats","text":"formats([io]; sortby=:format)\n\nDisplays in io (defaults to stdout if io is not given) a table with all formats supported by GeoIO.jl and the packages used to load and save each of them. \n\nOptionally, sort the table by the :extension, :load or :save columns using the sortby argument.\n\n\n\n\n\n","category":"function"},{"location":"data/","page":"Data","title":"Data","text":"using GeoIO\n\nzone = GeoIO.load(\"data/zone.shp\")","category":"page"},{"location":"data/#Artifacts","page":"Data","title":"Artifacts","text":"","category":"section"},{"location":"data/","page":"Data","title":"Data","text":"The GeoArtifacts.jl package provides utility functions to automatically download geospatial data from repositories on the internet.","category":"page"},{"location":"data/","page":"Data","title":"Data","text":"using GeoArtifacts\n\n# download artifacts from naturalearthdata.com\nearth = NaturalEarth.naturalearth1(\"water\")\nborders = NaturalEarth.borders()\nairports = NaturalEarth.airports()\nports = NaturalEarth.ports()\n\n# initialize viewer with a coarse \"raster\"\nearth |> Upscale(10, 5) |> viewer\n\n# add other elements to the visualization\nviz!(borders.geometry, color = \"cyan\")\nviz!(airports.geometry, color = \"black\", pointsize=4, pointmarker='✈')\nviz!(ports.geometry, color = \"blue\", pointsize=4, pointmarker='⚓')\n\n# display current figure\nMke.current_figure()","category":"page"},{"location":"visualization/#Visualization","page":"Visualization","title":"Visualization","text":"","category":"section"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"using GeoStats # hide\nimport CairoMakie as Mke # hide","category":"page"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"The framework provides powerful visualization recipes for geospatial data science via the Makie.jl project. These recipes were carefully designed to maximize productivity and to protect users from GIS jargon. The main entry point is the viz function:","category":"page"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"viz\nviz!","category":"page"},{"location":"visualization/#Meshes.viz","page":"Visualization","title":"Meshes.viz","text":"viz(object; [options])\n\nVisualize Meshes.jl object with various options.\n\nAvailable options\n\ncolor - color of geometries\nalpha - transparency in [0,1]\ncolormap - color scheme/map from ColorSchemes.jl\ncolorrange - minimum and maximum color values\nshowsegments - visualize segments\nsegmentcolor - color of segments\nsegmentsize - width of segments\nshowpoints - visualize points\npointmarker - marker of points\npointcolor - color of points\npointsize - size of points\n\nThe option color can be a single scalar or a vector of scalars. For Mesh subtypes, the length of the vector of colors determines if the colors should be assigned to vertices or to elements.\n\nExamples\n\nDifferent coloring methods (vertex vs. element):\n\n# vertex coloring (i.e. linear interpolation)\nviz(mesh, color = 1:nvertices(mesh))\n\n# element coloring (i.e. discrete colors)\nviz(mesh, color = 1:nelements(mesh))\n\nDifferent strategies to show the boundary of geometries (showsegments vs. boundary):\n\n# visualize boundary with showsegments\nviz(polygon, showsegments = true)\n\n# visualize boundary with separate call\nviz(polygon)\nviz!(boundary(polygon))\n\nNotes\n\nThis function will only work in the presence of a Makie.jl backend via package extensions in Julia v1.9 or later versions of the language.\n\n\n\n\n\n","category":"function"},{"location":"visualization/#Meshes.viz!","page":"Visualization","title":"Meshes.viz!","text":"viz!(object; [options])\n\nVisualize Meshes.jl object in an existing scene with options forwarded to viz.\n\n\n\n\n\n","category":"function"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"This function takes a geospatial domain as input and provides a set of aesthetic options to style the elements (i.e. geometries) of the domain.","category":"page"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"note: Note\nNotice that the geometry column of our geospatial data type is a domain (i.e. data.geometry isa Domain), and that this design enables several optimizations in the visualization itself.","category":"page"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"Users can also call Makie's plot function in the geometry column as in","category":"page"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"Mke.plot(data.geometry)","category":"page"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"and this is equivalent to calling the viz recipe above. The plot function also works with various other objects such as EmpiricalHistogram and EmpiricalVariogram. That is convenient if you don't remember the name of the recipe.","category":"page"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"Additionaly, we provide a basic scientific viewer to visualize all viewable variables in the data:","category":"page"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"viewer","category":"page"},{"location":"visualization/#GeoTables.viewer","page":"Visualization","title":"GeoTables.viewer","text":"viewer(geotable; kwargs...)\n\nBasic scientific viewer for geospatial table geotable.\n\nAesthetic options are forwarded via kwargs to the Meshes.viz recipe.\n\n\n\n\n\n","category":"function"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"Other plots are listed below that can be useful for geostatistical analysis.","category":"page"},{"location":"visualization/#Built-in","page":"Visualization","title":"Built-in","text":"","category":"section"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"A hscatter plot between two variables var1 and var2 (possibly with var2 = var1) is a simple scatter plot in which the dots represent all ordered pairs of values of var1 and var2 at a given lag h.","category":"page"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"𝒟 = georef((Z=[10sin(i/10) + j for i in 1:100, j in 1:200],))\n\n𝒮 = sample(𝒟, UniformSampling(500))\n\nfig = Mke.Figure(size = (800, 400))\nhscatter(fig[1,1], 𝒮, :Z, :Z, lag=0)\nhscatter(fig[1,2], 𝒮, :Z, :Z, lag=20)\nhscatter(fig[2,1], 𝒮, :Z, :Z, lag=40)\nhscatter(fig[2,2], 𝒮, :Z, :Z, lag=60)\nfig","category":"page"},{"location":"visualization/#PairPlots.jl","page":"Visualization","title":"PairPlots.jl","text":"","category":"section"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"The PairPlots.jl package provides the pairplot function that can be used with any table, including tables of attributes obtained with the values function.","category":"page"},{"location":"visualization/#Biplots.jl","page":"Visualization","title":"Biplots.jl","text":"","category":"section"},{"location":"visualization/","page":"Visualization","title":"Visualization","text":"The Biplot.jl package provides 2D and 3D statistical biplots.","category":"page"},{"location":"resources/education/#Education","page":"Education","title":"Education","text":"","category":"section"},{"location":"resources/education/","page":"Education","title":"Education","text":"We recommend the following educational resources.","category":"page"},{"location":"resources/education/#Learning-resources","page":"Education","title":"Learning resources","text":"","category":"section"},{"location":"resources/education/#Textbooks","page":"Education","title":"Textbooks","text":"","category":"section"},{"location":"resources/education/","page":"Education","title":"Education","text":"Hoffimann, J. 2023. Geospatial Data Science with Julia - A fresh approach to data science with geospatial data and the Julia programming language.\nIsaaks, E. and Srivastava, R. 1990. An Introduction to Applied Geostatistics - An introductory book on geostatistics that covers important concepts with very simple language.\nChilès, JP and Delfiner, P. 2012. Geostatistics: Modeling Spatial Uncertainty - An incredible book for those with good mathematical background.\nWebster, R and Oliver, MA. 2007. Geostatistics for Environmental Scientists - A great book with good balance between theory and practice.\nJournel, A. and Huijbregts, Ch. 2003. Mining Geostatistics - A great book with both theoretical and practical developments.","category":"page"},{"location":"resources/education/#Video-lectures","page":"Education","title":"Video lectures","text":"","category":"section"},{"location":"resources/education/","page":"Education","title":"Education","text":"Júlio Hoffimann - Video lectures with the GeoStats.jl framework.\nEdward Isaaks - Video lectures on variography, Kriging and related concepts.\nJef Caers - Video lectures on two-point and multiple-point methods.","category":"page"},{"location":"resources/education/#Workshop-material","page":"Education","title":"Workshop material","text":"","category":"section"},{"location":"resources/education/","page":"Education","title":"Education","text":"UFMG 2024 [Portuguese] - Geociência de Dados na Mineração, UFMG 2024\nJuliaEO 2024 [English] - Global Workshop on Earth Observation, AIRCentre 2024\nUFMG 2023 [Portuguese] - Geociência de Dados na Mineração, UFMG 2023\nJuliaEO 2023 [English] - Global Workshop on Earth Observation, AIRCentre 2023\nCBMina 2021 [Portuguese] - Introução à Geoestatística, CBMina 2021\nUFMG 2021 [Portuguese] - Introdução à Geoestatística, UFMG 2021","category":"page"},{"location":"resources/education/#Related-concepts","page":"Education","title":"Related concepts","text":"","category":"section"},{"location":"resources/education/#GaussianProcesses.jl","page":"Education","title":"GaussianProcesses.jl","text":"","category":"section"},{"location":"resources/education/","page":"Education","title":"Education","text":"GaussianProcesses.jl - Gaussian process regression and Simple Kriging are essentially two names for the same concept. The derivation of Kriging estimators, however; does not require distributional assumptions. It is a beautiful coincidence that for multivariate Gaussian distributions, Simple Kriging gives the conditional expectation.","category":"page"},{"location":"resources/education/#KernelFunctions.jl","page":"Education","title":"KernelFunctions.jl","text":"","category":"section"},{"location":"resources/education/","page":"Education","title":"Education","text":"KernelFunctions.jl - Spatial structure can be represented in many different forms: covariance, variogram, correlogram, etc. Variograms are more general than covariance kernels according to the intrinsic stationary property. This means that there are variogram models with no covariance counterpart. Furthermore, empirical variograms can be easily estimated from the data (in various directions) with an efficient procedure. GeoStats.jl treats variograms as first-class objects.","category":"page"},{"location":"resources/education/#Interpolations.jl","page":"Education","title":"Interpolations.jl","text":"","category":"section"},{"location":"resources/education/","page":"Education","title":"Education","text":"Interpolations.jl - Kriging and spline interpolation have different purposes, yet these two methods are sometimes listed as competing alternatives. Kriging estimation is about minimizing variance (or estimation error), whereas spline interpolation is about deriving smooth estimators for computer visualization. Kriging is a generalization of splines in which one has the freedom to customize geospatial structure based on data. Besides the estimate itself, Kriging also provides the variance map as a function of point patterns.","category":"page"},{"location":"resources/education/#ScikitLearn.jl","page":"Education","title":"ScikitLearn.jl","text":"","category":"section"},{"location":"resources/education/","page":"Education","title":"Education","text":"ScikitLearn.jl - Traditional statistical learning relies on core assumptions that do not hold in geospatial settings (fixed support, i.i.d. samples, ...). Geostatistical learning has been introduced recently as an attempt to push the frontiers of statistical learning with geospatial data.","category":"page"},{"location":"about/license/","page":"License","title":"License","text":"The GeoStats.jl project is licensed under the MIT license:","category":"page"},{"location":"about/license/","page":"License","title":"License","text":"MIT License\n\nCopyright (c) 2015 Júlio Hoffimann\n\n
","category":"page"},{"location":"transforms/#Geometric-transforms","page":"Transforms","title":"Geometric transforms","text":"","category":"section"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"Please check the Meshes.jl documentation for an updated list of geometric transforms. As an example consider the rotation of geospatial data over a Cartesian grid:","category":"page"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"# geospatial domain\nΩ = georef((Z=rand(10, 10),))\n\n# apply geometric transform\nΩr = Ω |> Rotate(Angle2d(π/4))\n\nfig = Mke.Figure(size = (800, 400))\nviz(fig[1,1], Ω.geometry, color = Ω.Z)\nviz(fig[1,2], Ωr.geometry, color = Ωr.Z)\nfig","category":"page"},{"location":"transforms/#Geostatistical-transforms","page":"Transforms","title":"Geostatistical transforms","text":"","category":"section"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"Bellow is the current list of transforms that operate on both the geometries and features of geospatial data. They are implemented in the GeoStatsBase.jl package.","category":"page"},{"location":"transforms/#UniqueCoords","page":"Transforms","title":"UniqueCoords","text":"","category":"section"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"UniqueCoords","category":"page"},{"location":"transforms/#GeoStatsTransforms.UniqueCoords","page":"Transforms","title":"GeoStatsTransforms.UniqueCoords","text":"UniqueCoords(var₁ => agg₁, var₂ => agg₂, ..., varₙ => aggₙ)\n\nRetain locations in data with unique coordinates.\n\nDuplicates of a variable varᵢ are aggregated with aggregation function aggᵢ. If an aggregation function is not defined for variable varᵢ, the default aggregation function will be used. Default aggregation function is mean for continuous variables and first otherwise.\n\nExamples\n\nUniqueCoords(1 => last, 2 => maximum)\nUniqueCoords(:a => first, :b => minimum)\nUniqueCoords(\"a\" => last, \"b\" => maximum)\n\n\n\n\n\n","category":"type"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"# point set with repeated points\np = rand(Point, 50)\nΩ = georef((Z=rand(100),), [p; p])","category":"page"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"# discard repeated points\n𝒰 = Ω |> UniqueCoords()","category":"page"},{"location":"transforms/#Detrend","page":"Transforms","title":"Detrend","text":"","category":"section"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"Detrend","category":"page"},{"location":"transforms/#GeoStatsTransforms.Detrend","page":"Transforms","title":"GeoStatsTransforms.Detrend","text":"Detrend(col₁, col₂, ..., colₙ; degree=1)\nDetrend([col₁, col₂, ..., colₙ]; degree=1)\nDetrend((col₁, col₂, ..., colₙ); degree=1)\n\nThe transform that detrends columns col₁, col₂, ..., colₙ with a polynomial of given degree.\n\nDetrend(regex; degree=1)\n\nDetrends the columns that match with regex.\n\nExamples\n\nDetrend(1, 3, 5)\nDetrend([:a, :c, :e])\nDetrend((\"a\", \"c\", \"e\"))\nDetrend(r\"[ace]\", degree=2)\nDetrend(:)\n\nReferences\n\nMenafoglio, A., Secchi, P. 2013. A Universal Kriging predictor for spatially dependent functional data of a Hilbert Space\n\n\n\n\n\n","category":"type"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"# quadratic trend + random noise\nr = range(-1, stop=1, length=100)\nμ = [x^2 + y^2 for x in r, y in r]\nϵ = 0.1rand(100, 100)\nΩ = georef((Z=μ+ϵ,))\n\n# detrend and obtain noise component\n𝒩 = Ω |> Detrend(:Z, degree=2)\n\nfig = Mke.Figure(size = (800, 400))\nviz(fig[1,1], Ω.geometry, color = Ω.Z)\nviz(fig[1,2], 𝒩.geometry, color = 𝒩.Z)\nfig","category":"page"},{"location":"transforms/#Potrace","page":"Transforms","title":"Potrace","text":"","category":"section"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"Potrace","category":"page"},{"location":"transforms/#GeoStatsTransforms.Potrace","page":"Transforms","title":"GeoStatsTransforms.Potrace","text":"Potrace(mask; [ϵ])\nPotrace(mask, var₁ => agg₁, ..., varₙ => aggₙ; [ϵ])\n\nTrace polygons on 2D image data with Selinger's Potrace algorithm.\n\nThe categories stored in column mask are converted into binary masks, which are then traced into multi-polygons. When provided, the option ϵ is forwarded to Selinger's simplification algorithm.\n\nDuplicates of a variable varᵢ are aggregated with aggregation function aggᵢ. If an aggregation function is not defined for variable varᵢ, the default aggregation function will be used. Default aggregation function is mean for continuous variables and first otherwise.\n\nExamples\n\nPotrace(:mask, ϵ=0.1)\nPotrace(1, 1 => last, 2 => maximum)\nPotrace(:mask, :a => first, :b => minimum)\nPotrace(\"mask\", \"a\" => last, \"b\" => maximum)\n\nReferences\n\nSelinger, P. 2003. Potrace: A polygon-based tracing algorithm\n\n\n\n\n\n","category":"type"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"# continuous feature\nZ = [sin(i/10) + sin(j/10) for i in 1:100, j in 1:100]\n\n# binary mask\nM = Z .> 0\n\n# georeference data\nΩ = georef((Z=Z, M=M))\n\n# trace polygons using mask\n𝒯 = Ω |> Potrace(:M)\n\nfig = Mke.Figure(size = (800, 400))\nviz(fig[1,1], Ω.geometry, color = Ω.Z)\nviz(fig[1,2], 𝒯.geometry, color = 𝒯.Z)\nfig","category":"page"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"𝒯.geometry","category":"page"},{"location":"transforms/#Rasterize","page":"Transforms","title":"Rasterize","text":"","category":"section"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"Rasterize","category":"page"},{"location":"transforms/#GeoStatsTransforms.Rasterize","page":"Transforms","title":"GeoStatsTransforms.Rasterize","text":"Rasterize(grid)\nRasterize(grid, var₁ => agg₁, ..., varₙ => aggₙ)\n\nRasterize geometries within specified grid.\n\nRasterize(nx, ny)\nRasterize(nx, ny, var₁ => agg₁, ..., varₙ => aggₙ)\n\nAlternatively, use the grid with size nx by ny obtained with discretization of the bounding box.\n\nDuplicates of a variable varᵢ are aggregated with aggregation function aggᵢ. If an aggregation function is not defined for variable varᵢ, the default aggregation function will be used. Default aggregation function is mean for continuous variables and first otherwise.\n\nExamples\n\ngrid = CartesianGrid(10, 10)\nRasterize(grid)\nRasterize(10, 10)\nRasterize(grid, 1 => last, 2 => maximum)\nRasterize(10, 10, 1 => last, 2 => maximum)\nRasterize(grid, :a => first, :b => minimum)\nRasterize(10, 10, :a => first, :b => minimum)\nRasterize(grid, \"a\" => last, \"b\" => maximum)\nRasterize(10, 10, \"a\" => last, \"b\" => maximum)\n\n\n\n\n\n","category":"type"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"A = [1, 2, 3, 4, 5]\nB = [1.1, 2.2, 3.3, 4.4, 5.5]\np1 = PolyArea((2, 0), (6, 2), (2, 2))\np2 = PolyArea((0, 6), (3, 8), (0, 10))\np3 = PolyArea((3, 6), (9, 6), (9, 9), (6, 9))\np4 = PolyArea((7, 0), (10, 0), (10, 4), (7, 4))\np5 = PolyArea((1, 3), (5, 3), (6, 6), (3, 8), (0, 6))\ngt = georef((; A, B), [p1, p2, p3, p4, p5])\n\nnt = gt |> Rasterize(20, 20)\n\nviz(nt.geometry, color = nt.A)","category":"page"},{"location":"transforms/#Clustering","page":"Transforms","title":"Clustering","text":"","category":"section"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"Unlike traditional clustering algorithms in machine learning, geostatistical clustering (a.k.a. domaining) algorithms consider both the features and the geospatial coordinates of the data.","category":"page"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"Consider the following data as an example:","category":"page"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"Ω = georef((Z=[10sin(i/10) + j for i in 1:4:100, j in 1:4:100],))\n\nviz(Ω.geometry, color = Ω.Z)","category":"page"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"GHC","category":"page"},{"location":"transforms/#GeoStatsTransforms.GHC","page":"Transforms","title":"GeoStatsTransforms.GHC","text":"GHC(k, λ; kern=:epanechnikov, link=:ward, as=:cluster)\n\nA transform for partitioning geospatial data into k clusters according to a range λ using Geostatistical Hierarchical Clustering (GHC). The larger the range the more connected are nearby samples.\n\nParameters\n\nk - Approximate number of clusters\nλ - Approximate range of kernel function in length units\nkern - Kernel function (:uniform, :triangular or :epanechnikov)\nlink - Linkage function (:single, :average, :complete, :ward or :ward_presquared)\nas - Variable name used to store clustering results\n\nReferences\n\nFouedjio, F. 2016. A hierarchical clustering method for multivariate geostatistical data\n\nNotes\n\nThe range parameter controls the sparsity pattern of the pairwise distances, which can greatly affect the computational performance of the GHC algorithm. We recommend choosing a range that is small enough to connect nearby samples. For example, clustering data over a 100x100 Cartesian grid with unit spacing is possible with λ=1.0 or λ=2.0 but the problem starts to become computationally unfeasible around λ=10.0 due to the density of points.\n\n\n\n\n\n","category":"type"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"𝒞 = Ω |> GHC(20, 1.0)\n\nviz(𝒞.geometry, color = 𝒞.cluster)","category":"page"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"GSC","category":"page"},{"location":"transforms/#GeoStatsTransforms.GSC","page":"Transforms","title":"GeoStatsTransforms.GSC","text":"GSC(k, m; σ=1.0, tol=1e-4, maxiter=10, weights=nothing, as=:cluster)\n\nA transform for partitioning geospatial data into k clusters using Geostatistical Spectral Clustering (GSC).\n\nParameters\n\nk - Desired number of clusters\nm - Multiplicative factor for adjacent weights\nσ - Standard deviation for exponential model (default to 1.0)\ntol - Tolerance of k-means algorithm (default to 1e-4)\nmaxiter - Maximum number of iterations (default to 10)\nweights - Dictionary with weights for each attribute (default to nothing)\nas - Variable name used to store clustering results\n\nReferences\n\nRomary et al. 2015. Unsupervised classification of multivariate geostatistical data: Two algorithms\n\nNotes\n\nThe algorithm implemented here is slightly different than the algorithm\n\ndescribed in Romary et al. 2015. Instead of setting Wᵢⱼ = 0 when i <-/-> j, we simply magnify the weight by a multiplicative factor Wᵢⱼ *= m when i <–> j. This leads to dense matrices but also better results in practice.\n\n\n\n\n\n","category":"type"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"𝒞 = Ω |> GSC(50, 2.0)\n\nviz(𝒞.geometry, color = 𝒞.cluster)","category":"page"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"SLIC","category":"page"},{"location":"transforms/#GeoStatsTransforms.SLIC","page":"Transforms","title":"GeoStatsTransforms.SLIC","text":"SLIC(k, m; tol=1e-4, maxiter=10, weights=nothing, as=:cluster)\n\nA transform for clustering geospatial data into approximately k clusters using Simple Linear Iterative Clustering (SLIC). The transform produces clusters of samples that are spatially connected based on a distance dₛ and that, at the same time, are similar in terms of vars with distance dᵥ. The tradeoff is controlled with a hyperparameter parameter m in an additive model dₜ = √(dᵥ² + m²(dₛ/s)²).\n\nParameters\n\nk - Approximate number of clusters\nm - Hyperparameter of SLIC model\ntol - Tolerance of k-means algorithm (default to 1e-4)\nmaxiter - Maximum number of iterations (default to 10)\nweights - Dictionary with weights for each attribute (default to nothing)\nas - Variable name used to store clustering results\n\nReferences\n\nAchanta et al. 2011. SLIC superpixels compared to state-of-the-art superpixel methods\n\n\n\n\n\n","category":"type"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"𝒞 = Ω |> SLIC(50, 0.01)\n\nviz(𝒞.geometry, color = 𝒞.cluster)","category":"page"},{"location":"transforms/#Interpolate","page":"Transforms","title":"Interpolate","text":"","category":"section"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"Interpolate\nInterpolateNeighbors\nInterpolateMissing\nInterpolateNaN","category":"page"},{"location":"transforms/#GeoStatsTransforms.Interpolate","page":"Transforms","title":"GeoStatsTransforms.Interpolate","text":"Interpolate(domain, vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])\nInterpolate([g₁, g₂, ..., gₙ], vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])\n\nInterpolate geospatial data on given domain or vector of geometries [g₁, g₂, ..., gₙ], using geostatistical models model₁, ..., modelₙ for variables vars₁, ..., varsₙ.\n\nInterpolate(domain, model=NN(); [parameters])\nInterpolate([g₁, g₂, ..., gₙ], model=NN(); [parameters])\n\nInterpolate geospatial data on given domain or vector of geometries [g₁, g₂, ..., gₙ], using geostatistical model for all variables.\n\nParameters\n\npoint - Perform interpolation on point support (default to true)\nprob - Perform probabilistic interpolation (default to false)\n\nSee also InterpolateNeighbors, InterpolateMissing, InterpolateNaN.\n\n\n\n\n\n","category":"type"},{"location":"transforms/#GeoStatsTransforms.InterpolateNeighbors","page":"Transforms","title":"GeoStatsTransforms.InterpolateNeighbors","text":"InterpolateNeighbors(domain, vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])\nInterpolateNeighbors([g₁, g₂, ..., gₙ], vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])\n\nInterpolate geospatial data on given domain or set of geometries g₁, g₂, ..., gₙ, using geostatistical models model₁, ..., modelₙ for variables vars₁, ..., varsₙ.\n\nInterpolateNeighbors(domain, model=NN(); [parameters])\nInterpolateNeighbors([g₁, g₂, ..., gₙ], model=NN(); [parameters])\n\nInterpolate geospatial data on given domain or set of geometries g₁, g₂, ..., gₙ, using geostatistical model for all variables.\n\nUnlike Interpolate, this transform uses neighbor search methods to fit geostatistical models at each interpolation location with a reduced number of measurements.\n\nParameters\n\nminneighbors - Minimum number of neighbors (default to 1)\nmaxneighbors - Maximum number of neighbors (default to 10)\nneighborhood - Search neighborhood (default to nothing)\ndistance - A distance defined in Distances.jl (default to Euclidean())\npoint - Perform interpolation on point support (default to true)\nprob - Perform probabilistic interpolation (default to false)\n\nThe maxneighbors parameter can be used to perform interpolation with a subset of measurements per prediction location. If maxneighbors is not provided, then all measurements are used.\n\nTwo neighborhood search methods are available:\n\nIf a neighborhood is provided, local prediction is performed by sliding the neighborhood in the domain.\nIf a neighborhood is not provided, the prediction is performed using maxneighbors nearest neighbors according to distance.\n\nSee also Interpolate, InterpolateMissing, InterpolateNaN.\n\n\n\n\n\n","category":"type"},{"location":"transforms/#GeoStatsTransforms.InterpolateMissing","page":"Transforms","title":"GeoStatsTransforms.InterpolateMissing","text":"InterpolateMissing(vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])\nInterpolateMissing(vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])\n\nInterpolate geospatial data on its own domain, using geostatistical models model₁, ..., modelₙ and non-missing values of the variables vars₁, ..., varsₙ.\n\nInterpolateMissing(model=NN(); [parameters])\nInterpolateMissing(model=NN(); [parameters])\n\nInterpolate geospatial data on its own domain, using geostatistical model and non-missing values of all variables.\n\nJust like InterpolateNeighbors, this transform uses neighbor search methods to fit geostatistical models at each interpolation location with a reduced number of measurements.\n\nParameters\n\nminneighbors - Minimum number of neighbors (default to 1)\nmaxneighbors - Maximum number of neighbors (default to 10)\nneighborhood - Search neighborhood (default to nothing)\ndistance - A distance defined in Distances.jl (default to Euclidean())\npoint - Perform interpolation on point support (default to true)\nprob - Perform probabilistic interpolation (default to false)\n\nThe maxneighbors parameter can be used to perform interpolation with a subset of measurements per prediction location. If maxneighbors is not provided, then all measurements are used.\n\nTwo neighborhood search methods are available:\n\nIf a neighborhood is provided, local prediction is performed by sliding the neighborhood in the domain.\nIf a neighborhood is not provided, the prediction is performed using maxneighbors nearest neighbors according to distance.\n\nSee also InterpolateNaN, InterpolateNeighbors, Interpolate.\n\n\n\n\n\n","category":"type"},{"location":"transforms/#GeoStatsTransforms.InterpolateNaN","page":"Transforms","title":"GeoStatsTransforms.InterpolateNaN","text":"InterpolateNaN(vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])\nInterpolateNaN(vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])\n\nInterpolate geospatial data on its own domain, using geostatistical models model₁, ..., modelₙ and non-NaN values of the variables vars₁, ..., varsₙ.\n\nInterpolateNaN(model=NN(); [parameters])\nInterpolateNaN(model=NN(); [parameters])\n\nInterpolate geospatial data on its own domain, using geostatistical model and non-NaN values of all variables.\n\nJust like InterpolateNeighbors, this transform uses neighbor search methods to fit geostatistical models at each interpolation location with a reduced number of measurements.\n\nParameters\n\nminneighbors - Minimum number of neighbors (default to 1)\nmaxneighbors - Maximum number of neighbors (default to 10)\nneighborhood - Search neighborhood (default to nothing)\ndistance - A distance defined in Distances.jl (default to Euclidean())\npoint - Perform interpolation on point support (default to true)\nprob - Perform probabilistic interpolation (default to false)\n\nThe maxneighbors parameter can be used to perform interpolation with a subset of measurements per prediction location. If maxneighbors is not provided, then all measurements are used.\n\nTwo neighborhood search methods are available:\n\nIf a neighborhood is provided, local prediction is performed by sliding the neighborhood in the domain.\nIf a neighborhood is not provided, the prediction is performed using maxneighbors nearest neighbors according to distance.\n\nSee also InterpolateMissing, InterpolateNeighbors, Interpolate.\n\n\n\n\n\n","category":"type"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"table = (; Z=[1.,0.,1.])\ncoord = [(25.,25.), (50.,75.), (75.,50.)]\ngeotable = georef(table, coord)\n\ngrid = CartesianGrid(100, 100)\n\nmodel = Kriging(GaussianVariogram(range=35.))\n\ninterp = geotable |> Interpolate(grid, model)\n\nviz(interp.geometry, color = interp.Z)","category":"page"},{"location":"transforms/#Simulate","page":"Transforms","title":"Simulate","text":"","category":"section"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"Simulate","category":"page"},{"location":"transforms/#GeoStatsTransforms.Simulate","page":"Transforms","title":"GeoStatsTransforms.Simulate","text":"Simulate(domain, vars₁ => process₁, ..., varsₙ => processₙ; [parameters])\nSimulate(domain, nreals, vars₁ => process₁, ..., varsₙ => processₙ; [parameters])\nSimulate([g₁, g₂, ..., gₙ], vars₁ => process₁, ..., varsₙ => processₙ; [parameters])\nSimulate([g₁, g₂, ..., gₙ], nreals, vars₁ => process₁, ..., varsₙ => processₙ; [parameters])\n\nSimulate nreals realizations of variables varsᵢ with geostatistical process processᵢ over given domain or vector of geometries [g₁, g₂, ..., gₙ].\n\nThe parameters are forwarded to the rand method of the geostatistical processes.\n\n\n\n\n\n","category":"type"},{"location":"transforms/#CookieCutter","page":"Transforms","title":"CookieCutter","text":"","category":"section"},{"location":"transforms/","page":"Transforms","title":"Transforms","text":"CookieCutter","category":"page"},{"location":"transforms/#GeoStatsTransforms.CookieCutter","page":"Transforms","title":"GeoStatsTransforms.CookieCutter","text":"CookieCutter(domain, parent => process, var₁ => procmap₁, ..., varₙ => procmapₙ; [parameters])\nCookieCutter(domain, nreals, parent => process, var₁ => procmap₁, ..., varₙ => procmapₙ; [parameters])\n\nSimulate nreals realizations of variable parent with geostatistical process process, and each child variable varsᵢ with process map procmapᵢ, over given domain.\n\nThe process map must be an iterable of pairs of the form: value => process. Each process in the map is related to a value of the parent realization, therefore the values of the child variables will be chosen according to the values of the corresponding parent realization.\n\nThe parameters are forwarded to the rand method of the geostatistical processes.\n\nExamples\n\nparent = QuiltingProcess(trainimg, (30, 30))\nchild0 = GaussianProcess(SphericalVariogram(range=20.0, sill=0.2))\nchild1 = GaussianProcess(SphericalVariogram(MetricBall((200.0, 20.0))))\ntransform = CookieCutter(domain, :parent => parent, :child => [0 => child0, 1 => child1])\n\n\n\n\n\n","category":"type"},{"location":"declustering/#Declustering","page":"Declustering","title":"Declustering","text":"","category":"section"},{"location":"declustering/","page":"Declustering","title":"Declustering","text":"Declustered statistics are statistics that make use of geospatial coordinates in an attempt to correct potential sampling bias:","category":"page"},{"location":"declustering/","page":"Declustering","title":"Declustering","text":"\n\n
","category":"page"},{"location":"declustering/","page":"Declustering","title":"Declustering","text":"The following statistics have geospatial semantics:","category":"page"},{"location":"declustering/","page":"Declustering","title":"Declustering","text":"mean(::AbstractGeoTable, ::Any)\nvar(::AbstractGeoTable, ::Any)\nquantile(::AbstractGeoTable, ::Any, ::Any)","category":"page"},{"location":"declustering/#Statistics.mean-Tuple{AbstractGeoTable, Any}","page":"Declustering","title":"Statistics.mean","text":"mean(data, v)\nmean(data, v, s)\n\nDeclustered mean of geospatial data. Optionally, specify the variable v and the block side s.\n\n\n\n\n\n","category":"method"},{"location":"declustering/#Statistics.var-Tuple{AbstractGeoTable, Any}","page":"Declustering","title":"Statistics.var","text":"var(data, v)\nvar(data, v, s)\n\nDeclustered variance of geospatial data. Optionally, specify the variable v and the block side s.\n\n\n\n\n\n","category":"method"},{"location":"declustering/#Statistics.quantile-Tuple{AbstractGeoTable, Any, Any}","page":"Declustering","title":"Statistics.quantile","text":"quantile(data, v, p)\nquantile(data, v, p, s)\n\nDeclustered quantile of geospatial data at probability p. Optionally, specify the variable v and the block side s.\n\n\n\n\n\n","category":"method"},{"location":"declustering/","page":"Declustering","title":"Declustering","text":"A histogram with geospatial semantics is also available where the heights of the bins are adjusted based on the coordinates of the samples:","category":"page"},{"location":"declustering/","page":"Declustering","title":"Declustering","text":"EmpiricalHistogram","category":"page"},{"location":"declustering/#GeoStatsBase.EmpiricalHistogram","page":"Declustering","title":"GeoStatsBase.EmpiricalHistogram","text":"EmpiricalHistogram(sdata, var, [s]; kwargs...)\n\nSpatial histogram of variable var in spatial data sdata. Optionally, specify the block side s and the keyword arguments kwargs for the fit(Histogram, ...) call.\n\n\n\n\n\n","category":"type"},{"location":"variography/fitting/#Fitting-variograms","page":"Fitting variograms","title":"Fitting variograms","text":"","category":"section"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"using GeoStats # hide\nimport CairoMakie as Mke # hide","category":"page"},{"location":"variography/fitting/#Overview","page":"Fitting variograms","title":"Overview","text":"","category":"section"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"Fitting theoretical variograms to empirical observations is an important modeling step to ensure valid mathematical models of geospatial continuity. Given an empirical variogram, the following function can be used to fit models:","category":"page"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"GeoStatsFunctions.fit","category":"page"},{"location":"variography/fitting/#GeoStatsFunctions.fit","page":"Fitting variograms","title":"GeoStatsFunctions.fit","text":"fit(F, f, algo=WeightedLeastSquares(); kwargs...)\n\nFit theoretical geostatistical function of type F to empirical function f using algorithm algo.\n\nOptionally fix theoretical parameters like range, sill and nugget in the kwargs.\n\nExamples\n\njulia> fit(SphericalVariogram, g)\njulia> fit(ExponentialVariogram, g)\njulia> fit(ExponentialVariogram, g, sill=1.0)\njulia> fit(ExponentialVariogram, g, maxsill=1.0)\njulia> fit(GaussianVariogram, g, WeightedLeastSquares())\n\n\n\n\n\nfit(Fs, f, algo=WeightedLeastSquares(); kwargs...)\n\nFit theoretical geostatistical functions of types Fs to empirical function f using algorithm algo and return the one with minimum error.\n\nExamples\n\njulia> fit([SphericalVariogram, ExponentialVariogram], g)\n\n\n\n\n\nfit(F, f, weightfun; kwargs...)\n\nConvenience method that forwards the weighting function weightfun to the WeightedLeastSquares algorithm.\n\nExamples\n\nfit(SphericalVariogram, g, h -> exp(-h))\nfit(Variogram, g, h -> exp(-h/100))\n\n\n\n\n\nfit(Variogram, g, algo=WeightedLeastSquares(); kwargs...)\n\nFit all stationary Variogram models to empirical variogram g using algorithm algo and return the one with minimum error.\n\nExamples\n\njulia> fit(Variogram, g)\njulia> fit(Variogram, g, WeightedLeastSquares())\n\n\n\n\n\n","category":"function"},{"location":"variography/fitting/#Example","page":"Fitting variograms","title":"Example","text":"","category":"section"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"# sinusoidal data\n𝒟 = georef((Z=[sin(i/2) + sin(j/2) for i in 1:50, j in 1:50],))\n\n# empirical variogram\ng = EmpiricalVariogram(𝒟, :Z, maxlag = 25u\"m\")\n\nvarioplot(g)","category":"page"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"We can fit specific models to the empirical variogram:","category":"page"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"γ = GeoStatsFunctions.fit(SineHoleVariogram, g)\n\nvarioplot(g)\nvarioplot!(γ, maxlag = 25u\"m\")\nMke.current_figure()","category":"page"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"or let the framework find the model with minimum error:","category":"page"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"γ = GeoStatsFunctions.fit(Variogram, g)\n\nvarioplot(g)\nvarioplot!(γ, maxlag = 25u\"m\")\nMke.current_figure()","category":"page"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"which should be a SineHoleVariogram given that the synthetic data of this example is sinusoidal.","category":"page"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"Optionally, we can specify a weighting function to give different weights to the lags:","category":"page"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"γ = GeoStatsFunctions.fit(SineHoleVariogram, g, h -> 1 / h^2)\n\nvarioplot(g)\nvarioplot!(γ, maxlag = 25u\"m\")\nMke.current_figure()","category":"page"},{"location":"variography/fitting/#Methods","page":"Fitting variograms","title":"Methods","text":"","category":"section"},{"location":"variography/fitting/#Weighted-least-squares","page":"Fitting variograms","title":"Weighted least squares","text":"","category":"section"},{"location":"variography/fitting/","page":"Fitting variograms","title":"Fitting variograms","text":"WeightedLeastSquares","category":"page"},{"location":"variography/fitting/#GeoStatsFunctions.WeightedLeastSquares","page":"Fitting variograms","title":"GeoStatsFunctions.WeightedLeastSquares","text":"WeightedLeastSquares()\nWeightedLeastSquares(w)\n\nFit theoretical variogram using weighted least squares with weighting function w (e.g. h -> 1/h). If no weighting function is provided, bin counts of empirical variogram are normalized and used as weights.\n\n\n\n\n\n","category":"type"},{"location":"quickstart/#Quickstart","page":"Quickstart","title":"Quickstart","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"A geostatistical workflow often consists of four steps:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Definition of geospatial data\nManipulation of geospatial data\nGeostatistical modeling\nScientific visualization","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"In this section, we walk through these steps to illustrate some of the features of the project. In the case of geostatistical modeling, we will specifically explore geostatistical learning models. If you prefer learning from video, check out the recording of our JuliaEO 2023 workshop:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"\n\n
","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"We assume that the following packages are loaded throughout the code examples:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"using GeoStats\nimport CairoMakie as Mke","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"The GeoStats.jl package exports the full stack for geospatial data science and geostatistical modeling. The CairoMakie.jl package is one of the possible visualization backends from the Makie.jl project.","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"If you are new to Julia and have never heard of Makie.jl before, here are a few tips to help you choose between the different backends:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"WGLMakie.jl is the preferred backend for interactive visualization on the web browser. It integrates well with Pluto.jl notebooks and other web-based applications.\nGLMakie.jl is the preferred backend for interactive high-performance visualization. It leverages native graphical resources and doesn't require a web browser to function.\nCairoMakie.jl is the preferred backend for publication-quality static visualization. It requires less computing power and is therefore recommended for those users with modest laptops.","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"note: Note\nWe recommend importing the Makie.jl backend as Mke to avoid polluting the Julia session with names from the visualization stack.","category":"page"},{"location":"quickstart/#Loading/creating-data","page":"Quickstart","title":"Loading/creating data","text":"","category":"section"},{"location":"quickstart/#Loading-data","page":"Quickstart","title":"Loading data","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"The Julia ecosystem for loading geospatial data is comprised of several low-level packages such as Shapefile.jl and GeoJSON.jl, which define their own very basic geometry types. Instead of requesting users to learn the so called GeoInterface.jl to handle these types, we provide the high-level GeoIO.jl package to load any file with geospatial data into well-tested geometries from the Meshes.jl submodule:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"using GeoIO\n\nzone = GeoIO.load(\"data/zone.shp\")","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Various functions are defined over these geometries, for instance:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"sum(area, zone.geometry)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"sum(perimeter, zone.geometry)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Please check the Meshes.jl documentation for more details.","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"note: Note\nWe highly recommend using Meshes.jl geometries in geospatial workflows as they were carefully designed to accomodate advanced features of the GeoStats.jl framework. Any other geometry type will likely fail with our geostatistical algorithms and pipelines.","category":"page"},{"location":"quickstart/#Creating-data","page":"Quickstart","title":"Creating data","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Geospatial data can be derived from other Julia variables. For example, given a Julia array, which is not attached to any coordinate system, we can georeference the array using the georef function:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Z = [10sin(i/10) + 2j for i in 1:50, j in 1:50]\n\nΩ = georef((Z=Z,))","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Default coordinates are assigned based on the size of the array, and different configurations can be obtained with different methods (see Data).","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Geospatial data can be visualized with the viz recipe function:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"viz(Ω.geometry, color = Ω.Z)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Alternatively, we provide a basic scientific viewer to visualize all viewable variables in the data with a colorbar and other interactive elements (see Visualization):","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"viewer(Ω)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"note: Note\nJulia supports unicode symbols with LaTeX syntax. We can type \\Omega and press TAB to get the symbol Ω in the example above. This autocompletion works in various text editors, including the VSCode editor with the Julia extension.","category":"page"},{"location":"quickstart/#Manipulating-data","page":"Quickstart","title":"Manipulating data","text":"","category":"section"},{"location":"quickstart/#Table-interface","page":"Quickstart","title":"Table interface","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Our geospatial data implements the Tables.jl interface, which means that they can be accessed as if they were tables with samples in the rows and variables in the columns. In this case, a special column named geometry is created on the fly, row by row, containing Meshes.jl geometries.","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"For those familiar with the productive DataFrames.jl interface, there is nothing new:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Ω[1,:]","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Ω[1,:geometry]","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Ω.Z","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"However, notice that our implementation performs some clever optimizations behind the scenes to avoid expensive creation of geometries:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Ω.geometry","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"We can always retrieve the table of attributes (or features) with the function values and the underlying geospatial domain with the function domain. This can be useful for writing algorithms that are type-stable and depend purely on the feature values:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"values(Ω)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"or on the geometries:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"domain(Ω)","category":"page"},{"location":"quickstart/#Geospatial-transforms","page":"Quickstart","title":"Geospatial transforms","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"It is easy to design advanced geospatial pipelines that operate on both the table of features and the underlying geospatial domain:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"# pipeline with table transforms\npipe = Quantile() → StdCoords()\n\n# feed geospatial data to pipeline\nΩ̂ = pipe(Ω)\n\n# plot distribution before and after pipeline\nfig = Mke.Figure(size = (800, 400))\nMke.hist(fig[1,1], Ω.Z, color = :gray)\nMke.hist(fig[2,1], Ω̂.Z, color = :gray)\nfig","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Coordinates before pipeline:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"boundingbox(Ω.geometry)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"and after pipeline:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"boundingbox(Ω̂.geometry)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"These pipelines are revertible meaning that we can transform the data, perform geostatistical modeling, and revert the pipelines to obtain estimates in the original sample space (see Transforms).","category":"page"},{"location":"quickstart/#Geospatial-queries","page":"Quickstart","title":"Geospatial queries","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"We provide three macros @groupby, @transform and @combine for powerful geospatial split-apply-combine patterns, as well as the function geojoin for advanced geospatial joins.","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"@transform(Ω, :W = 2 * :Z * area(:geometry))","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"These are very useful for geospatial data science as they hide the complexity of the geometry column. For more information, check the Queries section of the documentation.","category":"page"},{"location":"quickstart/#Geostatistical-modeling","page":"Quickstart","title":"Geostatistical modeling","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Having defined the geospatial data objects, we proceed and define the geostatistical learning model. Let's assume that we have geopatial data with some variable that we want to predict in a supervised learning setting. We load the data from a CSV file, and inspect the available columns:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"using CSV\nusing DataFrames\n\ntable = CSV.File(\"data/agriculture.csv\") |> DataFrame\n\nfirst(table, 5)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Columns band1, band2, ..., band4 represent four satellite bands for different locations (x, y) in this region. The column crop has the crop type for each location that was labeled manually with the purpose of fitting a learning model.","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"We can now georeference the table and plot some of the variables:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Ω = georef(table, (:x, :y))\n\nfig = Mke.Figure(size = (800, 400))\nviz(fig[1,1], Ω.geometry, color = Ω.band4)\nviz(fig[1,2], Ω.geometry, color = Ω.crop)\nfig","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Similar to a generic statistical learning workflow, we split the data into \"train\" and \"test\" sets. The main difference here is that our geospatial geosplit function accepts a separating plane specified by its normal direction (1,-1):","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Ωs, Ωt = geosplit(Ω, 0.2, (1.0, -1.0))\n\nviz(Ωs.geometry)\nviz!(Ωt.geometry, color = :gray90)\nMke.current_figure()","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"We can visualize the domain of the \"train\" (or source) set Ωs in blue, and the domain of the \"test\" (or target) set Ωt in gray. We reserved 20% of the samples to Ωs and 80% to Ωt. Internally, this geospatial geosplit function is implemented in terms of efficient geospatial partitions.","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Let's define our geostatistical learning model to predict the crop type based on the four satellite bands. We will use the DecisionTreeClassifier model, which is suitable for the task we want to perform. Any model from the StatsLeanModels.jl model is supported, including all models from ScikitLearn.jl:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"feats = [:band1, :band2, :band3, :band4]\nlabel = :crop\n\nmodel = DecisionTreeClassifier()","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"We will fit the model in Ωs where the features and labels are available and predict in Ωt where the features are available. The Learn transform automatically fits the model to the data:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"learn = Learn(Ωs, model, feats => label)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"The transform can be called with new data to generate predictions:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Ω̂t = learn(Ωt)","category":"page"},{"location":"quickstart/#Scientific-visualization","page":"Quickstart","title":"Scientific visualization","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"We note that the prediction of a geostatistical learning model is a geospatial data object, and we can inspect it with the same methods already described above. This also means that we can visualize the prediction directly, side by side with the true label in this synthetic example:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"fig = Mke.Figure(size = (800, 400))\nviz(fig[1,1], Ω̂t.geometry, color = Ω̂t.crop)\nviz(fig[1,2], Ωt.geometry, color = Ωt.crop)\nfig","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"With this example we conclude the basic workflow. To get familiar with other features of the project, please check the the reference guide.","category":"page"},{"location":"resources/ecosystem/#Ecosystem","page":"Ecosystem","title":"Ecosystem","text":"","category":"section"},{"location":"resources/ecosystem/","page":"Ecosystem","title":"Ecosystem","text":"The Julia ecosystem for geospatial modeling is maturing very quickly as the result of multiple initiatives such as JuliaEarth, JuliaClimate, and JuliaGeo. Each of these initiatives is associated with a different set of challenges that ultimatively determine the types of packages that are being developed in the corresponding GitHub organizations. In this section, we try to clarify what is available to first-time users of the language.","category":"page"},{"location":"resources/ecosystem/","page":"Ecosystem","title":"Ecosystem","text":"![](\"https://avatars.githubusercontent.com/u/59541313?s=200&v=4\")
\n![](\"https://avatars.githubusercontent.com/u/41747566?s=200&v=4\")
\n![](\"https://avatars.githubusercontent.com/u/10616454?s=200&v=4\")
","category":"page"},{"location":"resources/ecosystem/#JuliaEarth","page":"Ecosystem","title":"JuliaEarth","text":"","category":"section"},{"location":"resources/ecosystem/","page":"Ecosystem","title":"Ecosystem","text":"Originally created to host the GeoStats.jl framework, this initiative is primarily concerned with geospatial data science and geostatistical modeling. Due to the various applications in the subsurface of the Earth, most of our Julia packages were developed to work efficiently with both 2D and 3D geometries.","category":"page"},{"location":"resources/ecosystem/","page":"Ecosystem","title":"Ecosystem","text":"Unlike other initiatives, JuliaEarth is 100% Julia by design. This means that we do not rely on external libraries such as GDAL or Proj4 for geospatial work.","category":"page"},{"location":"resources/ecosystem/#JuliaClimate","page":"Ecosystem","title":"JuliaClimate","text":"","category":"section"},{"location":"resources/ecosystem/","page":"Ecosystem","title":"Ecosystem","text":"The most recent of the three initiatives, JuliaClimate has been created to address specific challenges in climate modeling. One of these challenges is access to climate data in Julia. Packages such as INMET.jl and CDSAPI.jl serve this purpose and are quite nice to work with.","category":"page"},{"location":"resources/ecosystem/#JuliaGeo","page":"Ecosystem","title":"JuliaGeo","text":"","category":"section"},{"location":"resources/ecosystem/","page":"Ecosystem","title":"Ecosystem","text":"Focused on bringing well-established external libraries to Julia, JuliaGeo provides packages that are widely used by geospatial communities from other programming languages. GDAL.jl, Proj4.jl and LibGEOS.jl are good examples of such packages.","category":"page"},{"location":"validation/#Validation","page":"Validation","title":"Validation","text":"","category":"section"},{"location":"validation/","page":"Validation","title":"Validation","text":"GeoStats.jl was designed to, among other things, facilitate rigorous comparison of different geostatistical models in the literature. As a user of geostatistics, you may be interested in trying various models on a given data set to pick the one with best performance. As a researcher in the field, you may be interested in benchmarking your new model against other established models.","category":"page"},{"location":"validation/","page":"Validation","title":"Validation","text":"Errors of geostatistical solvers can be estimated with the cverror function:","category":"page"},{"location":"validation/","page":"Validation","title":"Validation","text":"cverror","category":"page"},{"location":"validation/#GeoStatsValidation.cverror","page":"Validation","title":"GeoStatsValidation.cverror","text":"cverror(model::GeoStatsModel, geotable, method; kwargs...)\n\nEstimate error of model in a given geotable with error estimation method using Interpolate or InterpolateNeighbors depending on the passed kwargs.\n\ncverror(model::StatsLearnModel, geotable, method)\ncverror((model, invars => outvars), geotable, method)\n\nEstimate error of model in a given geotable with error estimation method using the Learn transform.\n\n\n\n\n\n","category":"function"},{"location":"validation/","page":"Validation","title":"Validation","text":"For example, we can perform block cross-validation on a decision tree model using the following code:","category":"page"},{"location":"validation/","page":"Validation","title":"Validation","text":"using GeoStats\nusing GeoIO\n\n# load geospatial data\nΩ = GeoIO.load(\"data/agriculture.csv\", coords = (\"x\", \"y\"))\n\n# 20%/80% split along the (1, -1) direction\nΩₛ, Ωₜ = geosplit(Ω, 0.2, (1.0, -1.0))\n\n# features and label for supervised learning\nfeats = [:band1,:band2,:band3,:band4]\nlabel = :crop\n\n# learning model\nmodel = DecisionTreeClassifier()\n\n# loss function\nloss = MisclassLoss()\n\n# block cross-validation with r = 30.\nbcv = BlockValidation(30., loss = Dict(:crop => loss))\n\n# estimate of generalization error\nϵ̂ = cverror((model, feats => label), Ωₛ, bcv)","category":"page"},{"location":"validation/","page":"Validation","title":"Validation","text":"We can unhide the labels in the target domain and compute the actual error for comparison:","category":"page"},{"location":"validation/","page":"Validation","title":"Validation","text":"# train in Ωₛ and predict in Ωₜ\nΩ̂ₜ = Ωₜ |> Learn(Ωₛ, model, feats => label)\n\t\n# actual error of the model\nϵ = mean(loss.(Ωₜ.crop, Ω̂ₜ.crop))","category":"page"},{"location":"validation/","page":"Validation","title":"Validation","text":"Below is the list of currently implemented validation methods.","category":"page"},{"location":"validation/#Leave-one-out","page":"Validation","title":"Leave-one-out","text":"","category":"section"},{"location":"validation/","page":"Validation","title":"Validation","text":"LeaveOneOut","category":"page"},{"location":"validation/#GeoStatsValidation.LeaveOneOut","page":"Validation","title":"GeoStatsValidation.LeaveOneOut","text":"LeaveOneOut(; loss=Dict())\n\nLeave-one-out validation. Optionally, specify loss function from LossFunctions.jl for some of the variables.\n\nReferences\n\nStone. 1974. Cross-Validatory Choice and Assessment of Statistical Predictions\n\n\n\n\n\n","category":"type"},{"location":"validation/#Leave-ball-out","page":"Validation","title":"Leave-ball-out","text":"","category":"section"},{"location":"validation/","page":"Validation","title":"Validation","text":"LeaveBallOut","category":"page"},{"location":"validation/#GeoStatsValidation.LeaveBallOut","page":"Validation","title":"GeoStatsValidation.LeaveBallOut","text":"LeaveBallOut(ball; loss=Dict())\n\nLeave-ball-out (a.k.a. spatial leave-one-out) validation. Optionally, specify loss function from the LossFunctions.jl package for some of the variables.\n\nLeaveBallOut(radius; loss=Dict())\n\nBy default, use Euclidean ball of given radius in space.\n\nReferences\n\nLe Rest et al. 2014. Spatial leave-one-out cross-validation for variable selection in the presence of spatial autocorrelation\n\n\n\n\n\n","category":"type"},{"location":"validation/#K-fold","page":"Validation","title":"K-fold","text":"","category":"section"},{"location":"validation/","page":"Validation","title":"Validation","text":"KFoldValidation","category":"page"},{"location":"validation/#GeoStatsValidation.KFoldValidation","page":"Validation","title":"GeoStatsValidation.KFoldValidation","text":"KFoldValidation(k; shuffle=true, loss=Dict())\n\nk-fold cross-validation. Optionally, shuffle the data, and specify loss function from LossFunctions.jl for some of the variables.\n\nReferences\n\nGeisser, S. 1975. The predictive sample reuse method with applications\nBurman, P. 1989. A comparative study of ordinary cross-validation, v-fold cross-validation and the repeated learning-testing methods\n\n\n\n\n\n","category":"type"},{"location":"validation/#Block","page":"Validation","title":"Block","text":"","category":"section"},{"location":"validation/","page":"Validation","title":"Validation","text":"BlockValidation","category":"page"},{"location":"validation/#GeoStatsValidation.BlockValidation","page":"Validation","title":"GeoStatsValidation.BlockValidation","text":"BlockValidation(sides; loss=Dict())\n\nCross-validation with blocks of given sides. Optionally, specify loss function from LossFunctions.jl for some of the variables. If only one side is provided, then blocks become cubes.\n\nReferences\n\nRoberts et al. 2017. Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure\nPohjankukka et al. 2017. Estimating the prediction performance of spatial models via spatial k-fold cross-validation\n\n\n\n\n\n","category":"type"},{"location":"validation/#Weighted","page":"Validation","title":"Weighted","text":"","category":"section"},{"location":"validation/","page":"Validation","title":"Validation","text":"WeightedValidation","category":"page"},{"location":"validation/#GeoStatsValidation.WeightedValidation","page":"Validation","title":"GeoStatsValidation.WeightedValidation","text":"WeightedValidation(weighting, folding; lambda=1.0, loss=Dict())\n\nAn error estimation method which samples are weighted with weighting method and split into folds with folding method. Weights are raised to lambda power in [0,1]. Optionally, specify loss function from LossFunctions.jl for some of the variables.\n\nReferences\n\nSugiyama et al. 2006. Importance-weighted cross-validation for covariate shift\nSugiyama et al. 2007. Covariate shift adaptation by importance weighted cross validation\n\n\n\n\n\n","category":"type"},{"location":"validation/#Density-ratio","page":"Validation","title":"Density-ratio","text":"","category":"section"},{"location":"validation/","page":"Validation","title":"Validation","text":"DensityRatioValidation","category":"page"},{"location":"validation/#GeoStatsValidation.DensityRatioValidation","page":"Validation","title":"GeoStatsValidation.DensityRatioValidation","text":"DensityRatioValidation(k; [parameters])\n\nDensity ratio validation where weights are first obtained with density ratio estimation, and then used in k-fold weighted cross-validation.\n\nParameters\n\nshuffle - Shuffle the data before folding (default to true)\nestimator - Density ratio estimator (default to LSIF())\noptlib - Optimization library (default to default_optlib(estimator))\nlambda - Power of density ratios (default to 1.0)\n\nPlease see DensityRatioEstimation.jl for a list of supported estimators.\n\nReferences\n\nHoffimann et al. 2020. Geostatistical Learning: Challenges and Opportunities\n\n\n\n\n\n","category":"type"},{"location":"queries/#Geospatial-queries","page":"Queries","title":"Geospatial queries","text":"","category":"section"},{"location":"queries/#Split-apply-combine","page":"Queries","title":"Split-apply-combine","text":"","category":"section"},{"location":"queries/","page":"Queries","title":"Queries","text":"We provide a geospatial version of the split-apply-combine pattern:","category":"page"},{"location":"queries/","page":"Queries","title":"Queries","text":"\n\n
","category":"page"},{"location":"queries/","page":"Queries","title":"Queries","text":"@groupby\n@transform\n@combine","category":"page"},{"location":"queries/#GeoTables.@groupby","page":"Queries","title":"GeoTables.@groupby","text":"@groupby(geotable, col₁, col₂, ..., colₙ)\n@groupby(geotable, [col₁, col₂, ..., colₙ])\n@groupby(geotable, (col₁, col₂, ..., colₙ))\n\nGroup geospatial geotable by columns col₁, col₂, ..., colₙ.\n\n@groupby(geotable, regex)\n\nGroup geospatial geotable by columns that match with regex.\n\nExamples\n\n@groupby(geotable, 1, 3, 5)\n@groupby(geotable, [:a, :c, :e])\n@groupby(geotable, (\"a\", \"c\", \"e\"))\n@groupby(geotable, r\"[ace]\")\n\n\n\n\n\n","category":"macro"},{"location":"queries/#GeoTables.@transform","page":"Queries","title":"GeoTables.@transform","text":"@transform(geotable, :col₁ = expr₁, :col₂ = expr₂, ..., :colₙ = exprₙ)\n\nReturns geospatial geotable with columns col₁, col₂, ..., colₙ computed with expressions expr₁, expr₂, ..., exprₙ.\n\nSee also: @groupby.\n\nExamples\n\n@transform(geotable, :z = :x + 2*:y)\n@transform(geotable, :w = :x^2 - :y^2)\n@transform(geotable, :sinx = sin(:x), :cosy = cos(:y))\n\ngroups = @groupby(geotable, :y)\n@transform(groups, :logx = log(:x))\n@transform(groups, :expz = exp(:z))\n\n@transform(geotable, {\"z\"} = {\"x\"} - 2*{\"y\"})\nxnm, ynm, znm = :x, :y, :z\n@transform(geotable, {znm} = {xnm} - 2*{ynm})\n\n\n\n\n\n","category":"macro"},{"location":"queries/#GeoTables.@combine","page":"Queries","title":"GeoTables.@combine","text":"@combine(geotable, :col₁ = expr₁, :col₂ = expr₂, ..., :colₙ = exprₙ)\n\nReturns geospatial geotable with columns :col₁, :col₂, ..., :colₙ computed with reduction expressions expr₁, expr₂, ..., exprₙ.\n\nIf a reduction expression is not defined for the :geometry column, the geometries will be reduced using Multi.\n\nSee also: @groupby.\n\nExamples\n\n@combine(geotable, :x_sum = sum(:x))\n@combine(geotable, :x_mean = mean(:x))\n@combine(geotable, :x_mean = mean(:x), :geometry = centroid(:geometry))\n\ngroups = @groupby(geotable, :y)\n@combine(groups, :x_prod = prod(:x))\n@combine(groups, :x_median = median(:x))\n@combine(groups, :x_median = median(:x), :geometry = centroid(:geometry))\n\n@combine(geotable, {\"z\"} = sum({\"x\"}) + prod({\"y\"}))\nxnm, ynm, znm = :x, :y, :z\n@combine(geotable, {znm} = sum({xnm}) + prod({ynm}))\n\n\n\n\n\n","category":"macro"},{"location":"queries/#Geospatial-join","page":"Queries","title":"Geospatial join","text":"","category":"section"},{"location":"queries/","page":"Queries","title":"Queries","text":"geojoin\ntablejoin","category":"page"},{"location":"queries/#GeoTables.geojoin","page":"Queries","title":"GeoTables.geojoin","text":"geojoin(geotable₁, geotable₂, var₁ => agg₁, ..., varₙ => aggₙ; kind=:left, pred=intersects, on=nothing)\n\nJoins geotable₁ with geotable₂ using a certain kind of join and predicate function pred that takes two geometries and returns a boolean ((g1, g2) -> g1 ⊆ g2).\n\nOptionally, add a variable value match to join, in addition to geometric match, by passing a single name or list of variable names (strings or symbols) to the on keyword argument. The variable value match use the isequal function.\n\nWhenever two or more matches are encountered, aggregate varᵢ with aggregation function aggᵢ. If no aggregation function is provided for a variable, then the aggregation function will be selected according to the scientific types: mean for continuous and first otherwise.\n\nKinds\n\n:left - Returns all rows of geotable₁ filling entries with missing when there is no match in geotable₂.\n:inner - Returns the subset of rows of geotable₁ that has a match in geotable₂.\n\nExamples\n\ngeojoin(gtb1, gtb2)\ngeojoin(gtb1, gtb2, 1 => mean)\ngeojoin(gtb1, gtb2, :a => mean, :b => std)\ngeojoin(gtb1, gtb2, \"a\" => mean, pred=issubset)\ngeojoin(gtb1, gtb2, on=:a)\ngeojoin(gtb1, gtb2, kind=:inner, on=[\"a\", \"b\"])\n\nSee also tablejoin.\n\n\n\n\n\n","category":"function"},{"location":"queries/#GeoTables.tablejoin","page":"Queries","title":"GeoTables.tablejoin","text":"tablejoin(geotable, table, var₁ => agg₁, ..., varₙ => aggₙ; kind=:left, on)\n\nJoins geotable with table using the values of the on variables to match rows with a certain kind of join.\n\non variables can be a single name or list of variable names (strings or symbols). The variable value match use the isequal function.\n\nWhenever two or more matches are encountered, aggregate varᵢ with aggregation function aggᵢ. If no aggregation function is provided for a variable, then the aggregation function will be selected according to the scientific types: mean for continuous and first otherwise.\n\nKinds\n\n:left - Returns all rows of geotable filling entries with missing when there is no match in table.\n:inner - Returns the subset of rows of geotable that has a match in table.\n\nExamples\n\ntablejoin(gtb, tab, on=:a)\ntablejoin(gtb, tab, 1 => mean, on=:a)\ntablejoin(gtb, tab, :a => mean, :b => std, on=:a)\ntablejoin(gtb, tab, \"a\" => mean, on=[:a, :b])\ntablejoin(gtb, tab, kind=:inner, on=[\"a\", \"b\"])\n\nSee also geojoin.\n\n\n\n\n\n","category":"function"},{"location":"","page":"Home","title":"Home","text":"GeoStats","category":"page"},{"location":"#GeoStats","page":"Home","title":"GeoStats","text":"(Image: GeoStats.jl)\n\n(Image: ) (Image: ) (Image: ) (Image: ) (Image: )\n\n(Image: ) (Image: )\n\nGeoStats.jl is an extensible framework for geospatial data science and geostatistical modeling fully written in Julia. It is comprised of several modules for advanced geometric processing, state-of-the-art geostatistical algorithms and sophisticated visualization of geospatial data.\n\nAll further information is provided in the online documentation.\n\n\n\n\n\n","category":"module"},{"location":"","page":"Home","title":"Home","text":"note: Star us on GitHub!\nIf you have found this software useful, please consider starring it on GitHub. This gives us an accurate lower bound of the (satisfied) user count.","category":"page"},{"location":"","page":"Home","title":"Home","text":"Organizations using the framework:","category":"page"},{"location":"","page":"Home","title":"Home","text":"\n ![](\"images/Petrobras.gif\")
\n ![](\"images/Vale.png\")
\n ![](\"images/Gazprom.png\")
\n ![](\"images/Nexa.jpg\")
\n ![](\"images/ENI.png\")
\n ![](\"images/Stanford.png\")
\n ![](\"images/UFMG.jpg\")
\n
\n \n ![](\"images/Arpeggeo.png\")
\n \n
\n\n
","category":"page"},{"location":"","page":"Home","title":"Home","text":"Consider reading the Geospatial Data Science with Julia book before reading this documentation. If you have questions, or would like to brainstorm ideas, don't hesitate to start a topic in our community channel.","category":"page"},{"location":"","page":"Home","title":"Home","text":"\n \n ![](\"https://juliaearth.github.io/geospatial-data-science-with-julia/images/cover.svg\")
\n \n
We provide a very powerful list of transforms that were designed to work seamlessly with geospatial data. These transforms are implemented in different packages depending on how they interact with geometries.
The list of supported transforms is continuously growing. The following code can be used to print an updated list in any project environment:
# packages to print type tree
+using InteractiveUtils
+using AbstractTrees
+using TransformsBase
+
+# packages with transforms
+using GeoStats
+
+# define the tree of types
+AbstractTrees.children(T::Type) = subtypes(T)
+
+# print all currently available transforms
+AbstractTrees.print_tree(TransformsBase.Transform)Transform
+├─ GeometricTransform
+│ ├─ Bridge
+│ ├─ CoordinateTransform
+│ │ ├─ Affine
+│ │ ├─ LengthUnit
+│ │ ├─ Morphological
+│ │ ├─ Proj
+│ │ ├─ Rotate
+│ │ ├─ Scale
+│ │ ├─ Stretch
+│ │ └─ Translate
+│ ├─ LambdaMuSmoothing
+│ ├─ Repair
+│ ├─ Shadow
+│ ├─ Slice
+│ └─ StdCoords
+├─ Identity
+├─ TableTransform
+│ ├─ Aggregate
+│ ├─ CookieCutter
+│ ├─ Detrend
+│ ├─ Downscale
+│ ├─ FeatureTransform
+│ │ ├─ EigenAnalysis
+│ │ ├─ Remainder
+│ │ ├─ ColwiseFeatureTransform
+│ │ │ ├─ Center
+│ │ │ ├─ Coalesce
+│ │ │ ├─ LowHigh
+│ │ │ ├─ Quantile
+│ │ │ └─ ZScore
+│ │ └─ StatelessFeatureTransform
+│ │ ├─ AbsoluteUnits
+│ │ ├─ Assert
+│ │ ├─ Closure
+│ │ ├─ Coerce
+│ │ ├─ ColTable
+│ │ ├─ Compose
+│ │ ├─ DropConstant
+│ │ ├─ DropExtrema
+│ │ ├─ DropMissing
+│ │ ├─ DropNaN
+│ │ ├─ DropUnits
+│ │ ├─ Filter
+│ │ ├─ Functional
+│ │ ├─ Indicator
+│ │ ├─ Learn
+│ │ ├─ Levels
+│ │ ├─ Map
+│ │ ├─ OneHot
+│ │ ├─ ProjectionPursuit
+│ │ ├─ Reject
+│ │ ├─ Rename
+│ │ ├─ Replace
+│ │ ├─ RowTable
+│ │ ├─ Sample
+│ │ ├─ Satisfies
+│ │ ├─ Select
+│ │ ├─ Sort
+│ │ ├─ StdFeats
+│ │ ├─ StdNames
+│ │ ├─ LogRatio
+│ │ │ ├─ ALR
+│ │ │ ├─ CLR
+│ │ │ └─ ILR
+│ │ ├─ Unit
+│ │ └─ Unitify
+│ ├─ ClusteringTransform
+│ │ ├─ GHC
+│ │ ├─ GSC
+│ │ └─ SLIC
+│ ├─ Interpolate
+│ ├─ InterpolateMissing
+│ ├─ InterpolateNaN
+│ ├─ InterpolateNeighbors
+│ ├─ Potrace
+│ ├─ Rasterize
+│ ├─ Simulate
+│ ├─ ParallelTableTransform
+│ ├─ Transfer
+│ ├─ UniqueCoords
+│ └─ Upscale
+└─ SequentialTransformTransforms at the leaves of the tree above should have a docstring with more information on the available options. For example, the documentation of the Select transform is shown below:
TableTransforms.Select — TypeSelect(col₁, col₂, ..., colₙ)
+Select([col₁, col₂, ..., colₙ])
+Select((col₁, col₂, ..., colₙ))The transform that selects columns col₁, col₂, ..., colₙ.
Select(col₁ => newcol₁, col₂ => newcol₂, ..., colₙ => newcolₙ)Selects the columns col₁, col₂, ..., colₙ and rename them to newcol₁, newcol₂, ..., newcolₙ.
Select(regex)Selects the columns that match with regex.
Examples
Select(1, 3, 5)
+Select([:a, :c, :e])
+Select(("a", "c", "e"))
+Select(1 => :x, 3 => :y)
+Select(:a => :x, :b => :y)
+Select("a" => "x", "b" => "y")
+Select(r"[ace]")Transforms of type FeatureTransform operate on the attribute table, whereas transforms of type GeometricTransform operate on the underlying geospatial domain:
TableTransforms.FeatureTransform — TypeFeatureTransformA transform that operates on the columns of the table containing features, i.e., simple attributes such as numbers, strings, etc.
Meshes.GeometricTransform — TypeGeometricTransformA method to transform the geometry (e.g. coordinates) of objects. See https://en.wikipedia.org/wiki/Geometric_transformation.
Other transforms such as Detrend are defined in terms of both the geospatial domain and the attribute table. All transforms and pipelines implement the following functions:
TransformsBase.isrevertible — Functionisrevertible(transform)Tells whether or not the transform is revertible, i.e. supports a revert function. Defaults to false for new transform types.
Transforms can be revertible and yet don't be invertible. Invertibility is a mathematical concept, whereas revertibility is a computational concept.
See also isinvertible.
TransformsBase.isinvertible — Functionisinvertible(transform)Tells whether or not the transform is invertible, i.e. whether it implements the inverse function. Defaults to false for new transform types.
Transforms can be invertible in the mathematical sense, i.e., there exists a one-to-one mapping between input and output spaces.
See also inverse, isrevertible.
TransformsBase.inverse — FunctionTransformsBase.apply — Functionnewobject, cache = apply(transform, object)Apply transform on the object. Return the newobject and a cache to revert the transform later.
TransformsBase.revert — Functionobject = revert(transform, newobject, cache)Revert the transform on the newobject using the cache from the corresponding apply call and return the original object. Only defined when the transform isrevertible.
TransformsBase.reapply — Functionnewobject = reapply(transform, object, cache)Reapply the transform to (a possibly different) object using a cache that was created with a previous apply call. Fallback to apply without using the cache.
Please check the TableTransforms.jl documentation for an updated list of feature transforms. As an example consider the following features over a Cartesian grid and their statistics:
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(Ω))| Row | variable | mean | min | median | max | nmissing | eltype |
|---|---|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Float64 | Float64 | Int64 | DataType | |
| 1 | a | 0.485276 | 0.00120786 | 0.484707 | 0.999245 | 0 | Float64 |
| 2 | b | 0.0176609 | -2.95863 | 0.0136271 | 3.62072 | 0 | Float64 |
| 3 | c | 0.50863 | 7.30319e-5 | 0.511608 | 0.999331 | 0 | Float64 |
We can create a pipeline that transforms the features to their normal quantile (or scores):
pipe = Quantile()
+
+Ω̄, cache = apply(pipe, Ω)
+
+describe(values(Ω̄))| Row | variable | mean | min | median | max | nmissing | eltype |
|---|---|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Float64 | Float64 | Int64 | DataType | |
| 1 | a | 0.00309023 | -3.09023 | 0.00125332 | 3.09023 | 0 | Float64 |
| 2 | b | 0.00309023 | -3.09023 | 0.00125332 | 3.09023 | 0 | Float64 |
| 3 | c | 0.00309023 | -3.09023 | 0.00125332 | 3.09023 | 0 | Float64 |
We can then revert the transform given any new geospatial data in the transformed sample space:
Ωₒ = revert(pipe, Ω̄, cache)
+
+describe(values(Ωₒ))| Row | variable | mean | min | median | max | nmissing | eltype |
|---|---|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Float64 | Float64 | Int64 | DataType | |
| 1 | a | 0.485759 | 0.00226717 | 0.48481 | 0.998913 | 0 | Float64 |
| 2 | b | 0.0206166 | -2.77883 | 0.0140045 | 3.60016 | 0 | Float64 |
| 3 | c | 0.509138 | 0.00012305 | 0.512828 | 0.998951 | 0 | Float64 |
The Learn transform is another important transform from StatsLearnModels.jl:
StatsLearnModels.Learn — TypeLearn(train, model, incols => outcols)Fits the statistical learning model using the input columns, selected by incols, and the output columns, selected by outcols, from the train table.
The column selection can be a single column identifier (index or name), a collection of identifiers or a regular expression (regex).
Examples
Learn(train, model, [1, 2, 3] => "d")
+Learn(train, model, [:a, :b, :c] => :d)
+Learn(train, model, ["a", "b", "c"] => 4)
+Learn(train, model, [1, 2, 3] => [:d, :e])
+Learn(train, model, r"[abc]" => ["d", "e"])For more details, consider watching our JuliaCon2021 talk:
+ +
Please check the Meshes.jl documentation for an updated list of geometric transforms. As an example consider the rotation of geospatial data over a Cartesian grid:
# geospatial domain
+Ω = georef((Z=rand(10, 10),))
+
+# apply geometric transform
+Ωr = Ω |> Rotate(Angle2d(π/4))
+
+fig = Mke.Figure(size = (800, 400))
+viz(fig[1,1], Ω.geometry, color = Ω.Z)
+viz(fig[1,2], Ωr.geometry, color = Ωr.Z)
+fig
Bellow is the current list of transforms that operate on both the geometries and features of geospatial data. They are implemented in the GeoStatsBase.jl package.
GeoStatsTransforms.UniqueCoords — TypeUniqueCoords(var₁ => agg₁, var₂ => agg₂, ..., varₙ => aggₙ)Retain locations in data with unique coordinates.
Duplicates of a variable varᵢ are aggregated with aggregation function aggᵢ. If an aggregation function is not defined for variable varᵢ, the default aggregation function will be used. Default aggregation function is mean for continuous variables and first otherwise.
Examples
UniqueCoords(1 => last, 2 => maximum)
+UniqueCoords(:a => first, :b => minimum)
+UniqueCoords("a" => last, "b" => maximum)# point set with repeated points
+p = rand(Point, 50)
+Ω = georef((Z=rand(100),), [p; p])| Z | +geometry | +
|---|---|
| Continuous | +Point | +
| [NoUnits] | +🖈 Cartesian{NoDatum} | +
| 0.0178461 | +(x: 0.992114 m, y: 0.524092 m, z: 0.301521 m) | +
| 0.573037 | +(x: 0.870208 m, y: 0.311693 m, z: 0.490476 m) | +
| 0.429644 | +(x: 0.441657 m, y: 0.392859 m, z: 0.577252 m) | +
| 0.153999 | +(x: 0.874132 m, y: 0.367076 m, z: 0.634166 m) | +
| 0.242723 | +(x: 0.459546 m, y: 0.0214615 m, z: 0.508143 m) | +
| 0.423039 | +(x: 0.741281 m, y: 0.240406 m, z: 0.258798 m) | +
| 0.703679 | +(x: 0.602122 m, y: 0.33676 m, z: 0.282085 m) | +
| 0.554076 | +(x: 0.837868 m, y: 0.269682 m, z: 0.381218 m) | +
| 0.977117 | +(x: 0.533105 m, y: 0.943335 m, z: 0.384398 m) | +
| 0.828324 | +(x: 0.815223 m, y: 0.793588 m, z: 0.149007 m) | +
| ⋮ | +⋮ | +
# discard repeated points
+𝒰 = Ω |> UniqueCoords()| Z | +geometry | +
|---|---|
| Continuous | +Point | +
| [NoUnits] | +🖈 Cartesian{NoDatum} | +
| 0.105255 | +(x: 0.992114 m, y: 0.524092 m, z: 0.301521 m) | +
| 0.750587 | +(x: 0.870208 m, y: 0.311693 m, z: 0.490476 m) | +
| 0.459804 | +(x: 0.441657 m, y: 0.392859 m, z: 0.577252 m) | +
| 0.180438 | +(x: 0.874132 m, y: 0.367076 m, z: 0.634166 m) | +
| 0.146086 | +(x: 0.459546 m, y: 0.0214615 m, z: 0.508143 m) | +
| 0.577357 | +(x: 0.741281 m, y: 0.240406 m, z: 0.258798 m) | +
| 0.82928 | +(x: 0.602122 m, y: 0.33676 m, z: 0.282085 m) | +
| 0.632146 | +(x: 0.837868 m, y: 0.269682 m, z: 0.381218 m) | +
| 0.622199 | +(x: 0.533105 m, y: 0.943335 m, z: 0.384398 m) | +
| 0.491146 | +(x: 0.815223 m, y: 0.793588 m, z: 0.149007 m) | +
| ⋮ | +⋮ | +
GeoStatsTransforms.Detrend — TypeDetrend(col₁, col₂, ..., colₙ; degree=1)
+Detrend([col₁, col₂, ..., colₙ]; degree=1)
+Detrend((col₁, col₂, ..., colₙ); degree=1)The transform that detrends columns col₁, col₂, ..., colₙ with a polynomial of given degree.
Detrend(regex; degree=1)Detrends the columns that match with regex.
Examples
Detrend(1, 3, 5)
+Detrend([:a, :c, :e])
+Detrend(("a", "c", "e"))
+Detrend(r"[ace]", degree=2)
+Detrend(:)References
# quadratic trend + random noise
+r = range(-1, stop=1, length=100)
+μ = [x^2 + y^2 for x in r, y in r]
+ϵ = 0.1rand(100, 100)
+Ω = georef((Z=μ+ϵ,))
+
+# detrend and obtain noise component
+𝒩 = Ω |> Detrend(:Z, degree=2)
+
+fig = Mke.Figure(size = (800, 400))
+viz(fig[1,1], Ω.geometry, color = Ω.Z)
+viz(fig[1,2], 𝒩.geometry, color = 𝒩.Z)
+fig
GeoStatsTransforms.Potrace — TypePotrace(mask; [ϵ])
+Potrace(mask, var₁ => agg₁, ..., varₙ => aggₙ; [ϵ])Trace polygons on 2D image data with Selinger's Potrace algorithm.
The categories stored in column mask are converted into binary masks, which are then traced into multi-polygons. When provided, the option ϵ is forwarded to Selinger's simplification algorithm.
Duplicates of a variable varᵢ are aggregated with aggregation function aggᵢ. If an aggregation function is not defined for variable varᵢ, the default aggregation function will be used. Default aggregation function is mean for continuous variables and first otherwise.
Examples
Potrace(:mask, ϵ=0.1)
+Potrace(1, 1 => last, 2 => maximum)
+Potrace(:mask, :a => first, :b => minimum)
+Potrace("mask", "a" => last, "b" => maximum)References
# 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
𝒯.geometry2 GeometrySet
+├─ Multi(4×PolyArea)
+└─ Multi(4×PolyArea)GeoStatsTransforms.Rasterize — TypeRasterize(grid)
+Rasterize(grid, var₁ => agg₁, ..., varₙ => aggₙ)Rasterize geometries within specified grid.
Rasterize(nx, ny)
+Rasterize(nx, ny, var₁ => agg₁, ..., varₙ => aggₙ)Alternatively, use the grid with size nx by ny obtained with discretization of the bounding box.
Duplicates of a variable varᵢ are aggregated with aggregation function aggᵢ. If an aggregation function is not defined for variable varᵢ, the default aggregation function will be used. Default aggregation function is mean for continuous variables and first otherwise.
Examples
grid = CartesianGrid(10, 10)
+Rasterize(grid)
+Rasterize(10, 10)
+Rasterize(grid, 1 => last, 2 => maximum)
+Rasterize(10, 10, 1 => last, 2 => maximum)
+Rasterize(grid, :a => first, :b => minimum)
+Rasterize(10, 10, :a => first, :b => minimum)
+Rasterize(grid, "a" => last, "b" => maximum)
+Rasterize(10, 10, "a" => last, "b" => maximum)A = [1, 2, 3, 4, 5]
+B = [1.1, 2.2, 3.3, 4.4, 5.5]
+p1 = PolyArea((2, 0), (6, 2), (2, 2))
+p2 = PolyArea((0, 6), (3, 8), (0, 10))
+p3 = PolyArea((3, 6), (9, 6), (9, 9), (6, 9))
+p4 = PolyArea((7, 0), (10, 0), (10, 4), (7, 4))
+p5 = PolyArea((1, 3), (5, 3), (6, 6), (3, 8), (0, 6))
+gt = georef((; A, B), [p1, p2, p3, p4, p5])
+
+nt = gt |> Rasterize(20, 20)
+
+viz(nt.geometry, color = nt.A)
Unlike traditional clustering algorithms in machine learning, geostatistical clustering (a.k.a. domaining) algorithms consider both the features and the geospatial coordinates of the data.
Consider the following data as an example:
Ω = georef((Z=[10sin(i/10) + j for i in 1:4:100, j in 1:4:100],))
+
+viz(Ω.geometry, color = Ω.Z)
GeoStatsTransforms.GHC — TypeGHC(k, λ; kern=:epanechnikov, link=:ward, as=:cluster)A transform for partitioning geospatial data into k clusters according to a range λ using Geostatistical Hierarchical Clustering (GHC). The larger the range the more connected are nearby samples.
Parameters
k - Approximate number of clustersλ - Approximate range of kernel function in length unitskern - Kernel function (:uniform, :triangular or :epanechnikov)link - Linkage function (:single, :average, :complete, :ward or :ward_presquared)as - Variable name used to store clustering resultsReferences
Notes
λ=1.0 or λ=2.0 but the problem starts to become computationally unfeasible around λ=10.0 due to the density of points.𝒞 = Ω |> GHC(20, 1.0)
+
+viz(𝒞.geometry, color = 𝒞.cluster)
GeoStatsTransforms.GSC — TypeGSC(k, m; σ=1.0, tol=1e-4, maxiter=10, weights=nothing, as=:cluster)A transform for partitioning geospatial data into k clusters using Geostatistical Spectral Clustering (GSC).
Parameters
k - Desired number of clustersm - Multiplicative factor for adjacent weightsσ - Standard deviation for exponential model (default to 1.0)tol - Tolerance of k-means algorithm (default to 1e-4)maxiter - Maximum number of iterations (default to 10)weights - Dictionary with weights for each attribute (default to nothing)as - Variable name used to store clustering resultsReferences
Notes
described in Romary et al. 2015. Instead of setting Wᵢⱼ = 0 when i <-/-> j, we simply magnify the weight by a multiplicative factor Wᵢⱼ *= m when i <–> j. This leads to dense matrices but also better results in practice.
𝒞 = Ω |> GSC(50, 2.0)
+
+viz(𝒞.geometry, color = 𝒞.cluster)
GeoStatsTransforms.SLIC — TypeSLIC(k, m; tol=1e-4, maxiter=10, weights=nothing, as=:cluster)A transform for clustering geospatial data into approximately k clusters using Simple Linear Iterative Clustering (SLIC). The transform produces clusters of samples that are spatially connected based on a distance dₛ and that, at the same time, are similar in terms of vars with distance dᵥ. The tradeoff is controlled with a hyperparameter parameter m in an additive model dₜ = √(dᵥ² + m²(dₛ/s)²).
Parameters
k - Approximate number of clustersm - Hyperparameter of SLIC modeltol - Tolerance of k-means algorithm (default to 1e-4)maxiter - Maximum number of iterations (default to 10)weights - Dictionary with weights for each attribute (default to nothing)as - Variable name used to store clustering resultsReferences
𝒞 = Ω |> SLIC(50, 0.01)
+
+viz(𝒞.geometry, color = 𝒞.cluster)
GeoStatsTransforms.Interpolate — TypeInterpolate(domain, vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])
+Interpolate([g₁, g₂, ..., gₙ], vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])Interpolate geospatial data on given domain or vector of geometries [g₁, g₂, ..., gₙ], using geostatistical models model₁, ..., modelₙ for variables vars₁, ..., varsₙ.
Interpolate(domain, model=NN(); [parameters])
+Interpolate([g₁, g₂, ..., gₙ], model=NN(); [parameters])Interpolate geospatial data on given domain or vector of geometries [g₁, g₂, ..., gₙ], using geostatistical model for all variables.
Parameters
point - Perform interpolation on point support (default to true)prob - Perform probabilistic interpolation (default to false)See also InterpolateNeighbors, InterpolateMissing, InterpolateNaN.
GeoStatsTransforms.InterpolateNeighbors — TypeInterpolateNeighbors(domain, vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])
+InterpolateNeighbors([g₁, g₂, ..., gₙ], vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])Interpolate geospatial data on given domain or set of geometries g₁, g₂, ..., gₙ, using geostatistical models model₁, ..., modelₙ for variables vars₁, ..., varsₙ.
InterpolateNeighbors(domain, model=NN(); [parameters])
+InterpolateNeighbors([g₁, g₂, ..., gₙ], model=NN(); [parameters])Interpolate geospatial data on given domain or set of geometries g₁, g₂, ..., gₙ, using geostatistical model for all variables.
Unlike Interpolate, this transform uses neighbor search methods to fit geostatistical models at each interpolation location with a reduced number of measurements.
Parameters
minneighbors - Minimum number of neighbors (default to 1)maxneighbors - Maximum number of neighbors (default to 10)neighborhood - Search neighborhood (default to nothing)distance - A distance defined in Distances.jl (default to Euclidean())point - Perform interpolation on point support (default to true)prob - Perform probabilistic interpolation (default to false)The maxneighbors parameter can be used to perform interpolation with a subset of measurements per prediction location. If maxneighbors is not provided, then all measurements are used.
Two neighborhood search methods are available:
If a neighborhood is provided, local prediction is performed by sliding the neighborhood in the domain.
If a neighborhood is not provided, the prediction is performed using maxneighbors nearest neighbors according to distance.
See also Interpolate, InterpolateMissing, InterpolateNaN.
GeoStatsTransforms.InterpolateMissing — TypeInterpolateMissing(vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])
+InterpolateMissing(vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])Interpolate geospatial data on its own domain, using geostatistical models model₁, ..., modelₙ and non-missing values of the variables vars₁, ..., varsₙ.
InterpolateMissing(model=NN(); [parameters])
+InterpolateMissing(model=NN(); [parameters])Interpolate geospatial data on its own domain, using geostatistical model and non-missing values of all variables.
Just like InterpolateNeighbors, this transform uses neighbor search methods to fit geostatistical models at each interpolation location with a reduced number of measurements.
Parameters
minneighbors - Minimum number of neighbors (default to 1)maxneighbors - Maximum number of neighbors (default to 10)neighborhood - Search neighborhood (default to nothing)distance - A distance defined in Distances.jl (default to Euclidean())point - Perform interpolation on point support (default to true)prob - Perform probabilistic interpolation (default to false)The maxneighbors parameter can be used to perform interpolation with a subset of measurements per prediction location. If maxneighbors is not provided, then all measurements are used.
Two neighborhood search methods are available:
If a neighborhood is provided, local prediction is performed by sliding the neighborhood in the domain.
If a neighborhood is not provided, the prediction is performed using maxneighbors nearest neighbors according to distance.
See also InterpolateNaN, InterpolateNeighbors, Interpolate.
GeoStatsTransforms.InterpolateNaN — TypeInterpolateNaN(vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])
+InterpolateNaN(vars₁ => model₁, ..., varsₙ => modelₙ; [parameters])Interpolate geospatial data on its own domain, using geostatistical models model₁, ..., modelₙ and non-NaN values of the variables vars₁, ..., varsₙ.
InterpolateNaN(model=NN(); [parameters])
+InterpolateNaN(model=NN(); [parameters])Interpolate geospatial data on its own domain, using geostatistical model and non-NaN values of all variables.
Just like InterpolateNeighbors, this transform uses neighbor search methods to fit geostatistical models at each interpolation location with a reduced number of measurements.
Parameters
minneighbors - Minimum number of neighbors (default to 1)maxneighbors - Maximum number of neighbors (default to 10)neighborhood - Search neighborhood (default to nothing)distance - A distance defined in Distances.jl (default to Euclidean())point - Perform interpolation on point support (default to true)prob - Perform probabilistic interpolation (default to false)The maxneighbors parameter can be used to perform interpolation with a subset of measurements per prediction location. If maxneighbors is not provided, then all measurements are used.
Two neighborhood search methods are available:
If a neighborhood is provided, local prediction is performed by sliding the neighborhood in the domain.
If a neighborhood is not provided, the prediction is performed using maxneighbors nearest neighbors according to distance.
See also InterpolateMissing, InterpolateNeighbors, Interpolate.
table = (; Z=[1.,0.,1.])
+coord = [(25.,25.), (50.,75.), (75.,50.)]
+geotable = georef(table, coord)
+
+grid = CartesianGrid(100, 100)
+
+model = Kriging(GaussianVariogram(range=35.))
+
+interp = geotable |> Interpolate(grid, model)
+
+viz(interp.geometry, color = interp.Z)
GeoStatsTransforms.Simulate — TypeSimulate(domain, vars₁ => process₁, ..., varsₙ => processₙ; [parameters])
+Simulate(domain, nreals, vars₁ => process₁, ..., varsₙ => processₙ; [parameters])
+Simulate([g₁, g₂, ..., gₙ], vars₁ => process₁, ..., varsₙ => processₙ; [parameters])
+Simulate([g₁, g₂, ..., gₙ], nreals, vars₁ => process₁, ..., varsₙ => processₙ; [parameters])Simulate nreals realizations of variables varsᵢ with geostatistical process processᵢ over given domain or vector of geometries [g₁, g₂, ..., gₙ].
The parameters are forwarded to the rand method of the geostatistical processes.
GeoStatsTransforms.CookieCutter — TypeCookieCutter(domain, parent => process, var₁ => procmap₁, ..., varₙ => procmapₙ; [parameters])
+CookieCutter(domain, nreals, parent => process, var₁ => procmap₁, ..., varₙ => procmapₙ; [parameters])Simulate nreals realizations of variable parent with geostatistical process process, and each child variable varsᵢ with process map procmapᵢ, over given domain.
The process map must be an iterable of pairs of the form: value => process. Each process in the map is related to a value of the parent realization, therefore the values of the child variables will be chosen according to the values of the corresponding parent realization.
The parameters are forwarded to the rand method of the geostatistical processes.
Examples
parent = QuiltingProcess(trainimg, (30, 30))
+child0 = GaussianProcess(SphericalVariogram(range=20.0, sill=0.2))
+child1 = GaussianProcess(SphericalVariogram(MetricBall((200.0, 20.0))))
+transform = CookieCutter(domain, :parent => parent, :child => [0 => child0, 1 => child1])Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
Transiograms of categorical variables are matrix-valued functions $t_{ab}(h)$ that measure the transition probability from categorical value $a$ to categorical value $b$ at a given lag $h \in \R$. They are often used for the simulation of Markov processes in more than one dimension.
The Carle's estimator of the empirical transiogram is given by
\[\widehat{t_{ab}}(h) = \frac{\sum_{(i,j) \in N(h)} {I_a}_i \cdot {I_b}_j}{\sum_{(i,j) \in N(h)} {I_a}_i}\]
where $N(h) = \left\{(i,j) \mid ||\x_i - \x_j|| = h\right\}$ is the set of pairs of locations at a distance $h$, and where $I_a$ and $I_b$ are the indicator variables for categorical values (or levels) $a$ and $b$, respectively.
Transiograms can be plotted with the following options:
GeoStatsFunctions.transioplot — Functiontransioplot(t; [options])Plot the transiogram t with given options.
Common transiogram options:
color - color of transiogramsize - size of transiogrammaxlag - maximum lag of variogramlevels - categorical levelsNotes
GeoStatsFunctions.EmpiricalTransiogram — TypeEmpiricalTransiogram(data, var; [parameters])Computes the empirical (a.k.a. experimental) omnidirectional transiogram for categorical variable var stored in geospatial data.
Parameters
20)Euclidean distance):ball)Available algorithms:
:full - loop over all pairs of points in the data:ball - loop over all points inside maximum lag ballAll implemented algorithms produce the exact same result. The :ball algorithm is considerably faster when the maximum lag is much smaller than the bounding box of the domain of the data.
See also: DirectionalTransiogram, PlanarTransiogram.
References
Carle, S.F. & Fogg, G.E. 1996. Transition probability-based indicator geostatistics
Carle et al 1998. Conditional Simulation of Hydrofacies Architecture: A Transition Probability/Markov Approach
GeoStatsFunctions.DirectionalTransiogram — FunctionDirectionalTransiogram(direction, data, var; dtol=1e-6u"m", [parameters])Computes the empirical transiogram for the categorical variable var stored in geospatial data along a given direction with band tolerance dtol in length units.
Optionally, forward parameters for the underlying EmpiricalTransiogram.
GeoStatsFunctions.PlanarTransiogram — FunctionPlanarTransiogram(normal, data, var; ntol=1e-6u"m", [parameters])Computes the empirical transiogram for the categorical variable var stored in geospatial data along a plane perpendicular to a normal direction with plane tolerance ntol in length units.
Optionally, forward parameters for the underlying EmpiricalTransiogram.
Consider the following categorical image:
using GeoStatsImages
+
+img = geostatsimage("Gaussian30x10")
+
+Z = [z < 0 ? 1 : 2 for z in img.Z]
+
+cat = georef((; Z=Z), img.geometry)
+
+cat |> viewer
We can estimate the ominidirectional transiogram with
t = EmpiricalTransiogram(cat, :Z, maxlag = 50.)
+
+transioplot(t)
Transiograms estimated along all directions in a given plane of reference are called transioplanes.
GeoStatsFunctions.EmpiricalTransioplane — TypeEmpiricalTransioplane(data, var;
+ normal=Vec(0,0,1), nangs=50,
+ ptol=0.5u"m", dtol=0.5u"m",
+ [parameters])Given a normal direction, estimate the transiogram of variable var along all directions in the corresponding plane of variation.
Optionally, specify the tolerance ptol in length units for the plane partition, the tolerance dtol in length units for the direction partition, the number of angles nangs in the plane, and forward the parameters to the underlying EmpiricalTransiogram.
The transioplane is plotted on a polar axis for all lags and angles:
t = EmpiricalTransioplane(cat, :Z, maxlag = 50.)
+
+planeplot(t)
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
We provide various theoretical transiogram models from the literature, which can can be combined with ellipsoid distances to model geometric anisotropy.
GeoStatsFunctions.LinearTransiogram — TypeLinearTransiogram(; ranges, proportions)
+LinearTransiogram(ball; proportions)A linear transiogram with ranges, and proportions. Optionally, use a custom metric ball.
transioplot(LinearTransiogram())
GeoStatsFunctions.GaussianTransiogram — TypeGaussianTransiogram(; ranges, proportions)
+GaussianTransiogram(ball; proportions)A Gaussian transiogram with ranges, and proportions. Optionally, use a custom metric ball.
transioplot(GaussianTransiogram())
GeoStatsFunctions.SphericalTransiogram — TypeSphericalTransiogram(; ranges, proportions)
+SphericalTransiogram(ball; proportions)A spherical transiogram with ranges, and proportions. Optionally, use a custom metric ball.
transioplot(SphericalTransiogram())
GeoStatsFunctions.ExponentialTransiogram — TypeExponentialTransiogram(; ranges, proportions)
+ExponentialTransiogram(ball; proportions)A exponential transiogram with ranges, and proportions. Optionally, use a custom metric ball.
transioplot(ExponentialTransiogram())
GeoStatsFunctions.MatrixExponentialTransiogram — TypeMatrixExponentialTransiogram(rate)
+MatrixExponentialTransiogram(ball, rate)An exponential transiogram with transition rate matrix. Optionally, specify a metric ball to model anisotropy.
MatrixExponentialTransiogram(lengths, proportions)
+MatrixExponentialTransiogram(ball, lengths, proportions)Alternatively, build transition rate matrix from mean lengths and relative proportions.
References
Carle, S.F. & Fogg, G.E. 1996. Transition probability-based indicator geostatistics
Carle et al 1998. Conditional Simulation of Hydrofacies Architecture: A Transition Probability/Markov Approach
transioplot(MatrixExponentialTransiogram((1.0u"m", 1.0u"m"), (0.5, 0.5)))
GeoStatsFunctions.PiecewiseLinearTransiogram — TypePiecewiseLinearTransiogram(abscissas, ordinates)
+PiecewiseLinearTransiogram(ball, abscissas, ordinates)A piecewise-linear transiogram model with abscissas and matrix ordinates obtained from an EmpiricalTransiogram. Optionally, specify a metric ball to model anisotropy.
References
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
GeoStats.jl was designed to, among other things, facilitate rigorous comparison of different geostatistical models in the literature. As a user of geostatistics, you may be interested in trying various models on a given data set to pick the one with best performance. As a researcher in the field, you may be interested in benchmarking your new model against other established models.
Errors of geostatistical solvers can be estimated with the cverror function:
GeoStatsValidation.cverror — Functioncverror(model::GeoStatsModel, geotable, method; kwargs...)Estimate error of model in a given geotable with error estimation method using Interpolate or InterpolateNeighbors depending on the passed kwargs.
cverror(model::StatsLearnModel, geotable, method)
+cverror((model, invars => outvars), geotable, method)Estimate error of model in a given geotable with error estimation method using the Learn transform.
For example, we can perform block cross-validation on a decision tree model using the following code:
using GeoStats
+using GeoIO
+
+# load geospatial data
+Ω = GeoIO.load("data/agriculture.csv", coords = ("x", "y"))
+
+# 20%/80% split along the (1, -1) direction
+Ωₛ, Ωₜ = geosplit(Ω, 0.2, (1.0, -1.0))
+
+# features and label for supervised learning
+feats = [:band1,:band2,:band3,:band4]
+label = :crop
+
+# learning model
+model = DecisionTreeClassifier()
+
+# loss function
+loss = MisclassLoss()
+
+# block cross-validation with r = 30.
+bcv = BlockValidation(30., loss = Dict(:crop => loss))
+
+# estimate of generalization error
+ϵ̂ = cverror((model, feats => label), Ωₛ, bcv)Dict{Symbol, Float64} with 1 entry:
+ :crop => 0.231826We can unhide the labels in the target domain and compute the actual error for comparison:
# train in Ωₛ and predict in Ωₜ
+Ω̂ₜ = Ωₜ |> Learn(Ωₛ, model, feats => label)
+
+# actual error of the model
+ϵ = mean(loss.(Ωₜ.crop, Ω̂ₜ.crop))0.23606322557004614Below is the list of currently implemented validation methods.
GeoStatsValidation.LeaveOneOut — TypeLeaveOneOut(; loss=Dict())Leave-one-out validation. Optionally, specify loss function from LossFunctions.jl for some of the variables.
References
GeoStatsValidation.LeaveBallOut — TypeLeaveBallOut(ball; loss=Dict())Leave-ball-out (a.k.a. spatial leave-one-out) validation. Optionally, specify loss function from the LossFunctions.jl package for some of the variables.
LeaveBallOut(radius; loss=Dict())By default, use Euclidean ball of given radius in space.
References
GeoStatsValidation.KFoldValidation — TypeKFoldValidation(k; shuffle=true, loss=Dict())k-fold cross-validation. Optionally, shuffle the data, and specify loss function from LossFunctions.jl for some of the variables.
References
GeoStatsValidation.BlockValidation — TypeBlockValidation(sides; loss=Dict())Cross-validation with blocks of given sides. Optionally, specify loss function from LossFunctions.jl for some of the variables. If only one side is provided, then blocks become cubes.
References
GeoStatsValidation.WeightedValidation — TypeWeightedValidation(weighting, folding; lambda=1.0, loss=Dict())An error estimation method which samples are weighted with weighting method and split into folds with folding method. Weights are raised to lambda power in [0,1]. Optionally, specify loss function from LossFunctions.jl for some of the variables.
References
GeoStatsValidation.DensityRatioValidation — TypeDensityRatioValidation(k; [parameters])Density ratio validation where weights are first obtained with density ratio estimation, and then used in k-fold weighted cross-validation.
Parameters
shuffle - Shuffle the data before folding (default to true)estimator - Density ratio estimator (default to LSIF())optlib - Optimization library (default to default_optlib(estimator))lambda - Power of density ratios (default to 1.0)Please see DensityRatioEstimation.jl for a list of supported estimators.
References
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
Variograms are widely used in geostatistics due to their intimate connection with (cross-)variance and visual interpretability. The following video explains the concept in detail:
+ +
The Matheron's estimator of the empirical variogram is given by
\[\widehat{\gamma_M}(h) = \frac{1}{2|N(h)|} \sum_{(i,j) \in N(h)} (z_i - z_j)^2\]
where $N(h) = \left\{(i,j) \mid ||\x_i - \x_j|| = h\right\}$ is the set of pairs of locations at a distance $h$ and $|N(h)|$ is the cardinality of the set. Alternatively, the robust Cressie's estimator is given by
\[\widehat{\gamma_C}(h) = \frac{1}{2}\frac{\left\{\frac{1}{|N(h)|} \sum_{(i,j) \in N(h)} |z_i - z_j|^{1/2}\right\}^4}{0.457 + \frac{0.494}{|N(h)|} + \frac{0.045}{|N(h)|^2}}\]
Both estimators are available and can be used with general distance functions in order to for example:
Please see Distances.jl for a complete list of distance functions.
The high-performance estimation procedure implemented in the framework can consider all pairs of locations regardless of direction (ominidirectional) or a specified partition of the geospatial data (e.g. directional, planar).
Variograms can plotted with the following options:
GeoStatsFunctions.varioplot — Functionvarioplot(γ; [options])Plot the variogram γ with given options.
Common variogram options:
color - color of variogramsize - size of variogrammaxlag - maximum lag of variogramEmpirical variogram options:
pointsize - size of points of variogramshowtext - show text countstextsize - size of text countsshowhist - show histogramhistcolor - color of histogramNotes
GeoStatsFunctions.EmpiricalVariogram — TypeEmpiricalVariogram(data, var₁, var₂=var₁; [parameters])Computes the empirical (a.k.a. experimental) omnidirectional (cross-)variogram for variables var₁ and var₂ stored in geospatial data.
Parameters
20)Euclidean distance):matheron estimator):ball)Available estimators:
:matheron - simple estimator based on squared differences:cressie - robust estimator based on 4th power of differencesAvailable algorithms:
:full - loop over all pairs of points in the data:ball - loop over all points inside maximum lag ballAll implemented algorithms produce the exact same result. The :ball algorithm is considerably faster when the maximum lag is much smaller than the bounding box of the domain of the data.
See also: DirectionalVariogram, PlanarVariogram, EmpiricalVarioplane.
References
Chilès, JP and Delfiner, P. 2012. Geostatistics: Modeling Spatial Uncertainty
Webster, R and Oliver, MA. 2007. Geostatistics for Environmental Scientists
Hoffimann, J and Zadrozny, B. 2019. Efficient variography with partition variograms
GeoStatsFunctions.DirectionalVariogram — FunctionDirectionalVariogram(direction, data, var₁, var₂=var₁; dtol=1e-6u"m", [parameters])Computes the empirical (cross-)variogram for the variables var₁ and var₂ stored in geospatial data along a given direction with band tolerance dtol in length units.
Optionally, forward parameters for the underlying EmpiricalVariogram.
GeoStatsFunctions.PlanarVariogram — FunctionPlanarVariogram(normal, data, var₁, var₂=var₁; ntol=1e-6u"m", [parameters])Computes the empirical (cross-)variogram for the variables var₁ and var₂ stored in geospatial data along a plane perpendicular to a normal direction with plane tolerance ntol in length units.
Optionally, forward parameters for the underlying EmpiricalVariogram.
Consider the following example image:
using GeoStatsImages
+
+img = geostatsimage("Gaussian30x10")
+
+img |> viewer
We can estimate ominidirectional variograms, which consider pairs of points along all directions:
γ = EmpiricalVariogram(img, :Z, maxlag = 50.)
+
+varioplot(γ)
directional variograms along a specific direction:
γₕ = DirectionalVariogram((1.,0.), img, :Z, maxlag = 50.)
+γᵥ = DirectionalVariogram((0.,1.), img, :Z, maxlag = 50.)
+
+varioplot(γₕ, color = :maroon, histcolor = :maroon)
+varioplot!(γᵥ)
+Mke.current_figure()
or planar variograms over a specific plane:
γᵥ = PlanarVariogram((1.,0.), img, :Z, maxlag = 50.)
+γₕ = PlanarVariogram((0.,1.), img, :Z, maxlag = 50.)
+
+varioplot(γₕ, color = :maroon, histcolor = :maroon)
+varioplot!(γᵥ)
+Mke.current_figure()
The directional and planar variograms coincide in this example because planes are equal to lines in 2-dimensional space. These concepts are most useful in 3-dimensional space where we may be interested in comparing the horizontal planar range to the vertical directional range.
Variograms estimated along all directions in a given plane of reference are called varioplanes.
GeoStatsFunctions.EmpiricalVarioplane — TypeEmpiricalVarioplane(data, var₁, var₂=var₁;
+ normal=Vec(0,0,1), nangs=50,
+ ptol=0.5u"m", dtol=0.5u"m",
+ [parameters])Given a normal direction, estimate the (cross-)variogram of variables var₁ and var₂ along all directions in the corresponding plane of variation.
Optionally, specify the tolerance ptol in length units for the plane partition, the tolerance dtol in length units for the direction partition, the number of angles nangs in the plane, and forward the parameters to the underlying EmpiricalVariogram.
The varioplane is plotted on a polar axis for all lags and angles:
γ = EmpiricalVarioplane(img, :Z, maxlag = 50.)
+
+planeplot(γ)
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
Fitting theoretical variograms to empirical observations is an important modeling step to ensure valid mathematical models of geospatial continuity. Given an empirical variogram, the following function can be used to fit models:
GeoStatsFunctions.fit — Functionfit(F, f, algo=WeightedLeastSquares(); kwargs...)Fit theoretical geostatistical function of type F to empirical function f using algorithm algo.
Optionally fix theoretical parameters like range, sill and nugget in the kwargs.
Examples
julia> fit(SphericalVariogram, g)
+julia> fit(ExponentialVariogram, g)
+julia> fit(ExponentialVariogram, g, sill=1.0)
+julia> fit(ExponentialVariogram, g, maxsill=1.0)
+julia> fit(GaussianVariogram, g, WeightedLeastSquares())fit(Fs, f, algo=WeightedLeastSquares(); kwargs...)Fit theoretical geostatistical functions of types Fs to empirical function f using algorithm algo and return the one with minimum error.
Examples
julia> fit([SphericalVariogram, ExponentialVariogram], g)fit(F, f, weightfun; kwargs...)Convenience method that forwards the weighting function weightfun to the WeightedLeastSquares algorithm.
Examples
fit(SphericalVariogram, g, h -> exp(-h))
+fit(Variogram, g, h -> exp(-h/100))fit(Variogram, g, algo=WeightedLeastSquares(); kwargs...)Fit all stationary Variogram models to empirical variogram g using algorithm algo and return the one with minimum error.
Examples
julia> fit(Variogram, g)
+julia> fit(Variogram, g, WeightedLeastSquares())# sinusoidal data
+𝒟 = georef((Z=[sin(i/2) + sin(j/2) for i in 1:50, j in 1:50],))
+
+# empirical variogram
+g = EmpiricalVariogram(𝒟, :Z, maxlag = 25u"m")
+
+varioplot(g)
We can fit specific models to the empirical variogram:
γ = GeoStatsFunctions.fit(SineHoleVariogram, g)
+
+varioplot(g)
+varioplot!(γ, maxlag = 25u"m")
+Mke.current_figure()
or let the framework find the model with minimum error:
γ = GeoStatsFunctions.fit(Variogram, g)
+
+varioplot(g)
+varioplot!(γ, maxlag = 25u"m")
+Mke.current_figure()
which should be a SineHoleVariogram given that the synthetic data of this example is sinusoidal.
Optionally, we can specify a weighting function to give different weights to the lags:
γ = GeoStatsFunctions.fit(SineHoleVariogram, g, h -> 1 / h^2)
+
+varioplot(g)
+varioplot!(γ, maxlag = 25u"m")
+Mke.current_figure()
GeoStatsFunctions.WeightedLeastSquares — TypeWeightedLeastSquares()
+WeightedLeastSquares(w)Fit theoretical variogram using weighted least squares with weighting function w (e.g. h -> 1/h). If no weighting function is provided, bin counts of empirical variogram are normalized and used as weights.
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
We provide various theoretical variogram models from the literature, which can can be combined with ellipsoid distances to model geometric anisotropy and nested with scalars or matrix coefficients to express multivariate relations.
In an intrinsic isotropic model, the variogram is only a function of the distance between any two points $\x_1,\x_2 \in \R^m$:
\[\gamma(\x_1,\x_2) = \gamma(||\x_1 - \x_2||) = \gamma(h)\]
Under the additional assumption of 2nd-order stationarity, the well-known covariance is directly related via $\gamma(h) = \cov(0) - \cov(h)$. This package implements a few commonly used as well as other more exotic variogram models. Most of these models share a set of default parameters (e.g. sill=1.0, range=1.0), which can be set with keyword arguments.
Functions are provided to query properties of variogram models programmatically:
GeoStatsFunctions.isstationary — Methodisstationary(f)Check if geostatistical function f possesses the 2nd-order stationary property.
Meshes.isisotropic — Methodisisotropic(f)Tells whether or not the geostatistical function f is isotropic.
Corresponding covariance models are available for convenience. For example, the GaussianVariogram and the GaussianCovariance are two alternative representations of the same geostatistical behavior.
\[\gamma(h) = (s - n) \left[1 - \exp\left(-3\left(\frac{h}{r}\right)^2\right)\right] + n \cdot \1_{(0,\infty)}(h)\]
GeoStatsFunctions.GaussianVariogram — TypeGaussianVariogram(; range, sill, nugget)
+GaussianVariogram(ball; sill, nugget)A Gaussian variogram with range, sill and nugget. Optionally, use a custom metric ball.
varioplot(GaussianVariogram())
\[\gamma(h) = (s - n) \left[\left(\frac{3}{2}\left(\frac{h}{r}\right) + \frac{1}{2}\left(\frac{h}{r}\right)^3\right) \cdot \1_{(0,r)}(h) + \1_{[r,\infty)}(h)\right] + n \cdot \1_{(0,\infty)}(h)\]
GeoStatsFunctions.SphericalVariogram — TypeSphericalVariogram(; range, sill, nugget)
+SphericalVariogram(ball; sill, nugget)A spherical variogram with range, sill and nugget. Optionally, use a custom metric ball.
varioplot(SphericalVariogram())
\[\gamma(h) = (s - n) \left[1 - \exp\left(-3\left(\frac{h}{r}\right)\right)\right] + n \cdot \1_{(0,\infty)}(h)\]
GeoStatsFunctions.ExponentialVariogram — TypeExponentialVariogram(; range, sill, nugget)
+ExponentialVariogram(ball; sill, nugget)A exponential variogram with range, sill and nugget. Optionally, use a custom metric ball.
varioplot(ExponentialVariogram())
\[\gamma(h) = (s - n) \left[1 - \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\sqrt{2\nu}\ 3\left(\frac{h}{r}\right)\right)^\nu K_\nu\left(\sqrt{2\nu}\ 3\left(\frac{h}{r}\right)\right)\right] + n \cdot \1_{(0,\infty)}(h)\]
GeoStatsFunctions.MaternVariogram — TypeMaternVariogram(; range, sill, nugget, order)
+MaternVariogram(ball; sill, nugget, order)A Matérn variogram with range, sill and nugget. Optionally, use a custom metric ball and order of the Bessel function.
varioplot(MaternVariogram())
\[\gamma(h) = (s - n) \left[\left(7\left(\frac{h}{r}\right)^2 - \frac{35}{4}\left(\frac{h}{r}\right)^3 + \frac{7}{2}\left(\frac{h}{r}\right)^5 - \frac{3}{4}\left(\frac{h}{r}\right)^7\right) \cdot \1_{(0,r)}(h) + \1_{[r,\infty)}(h)\right] + n \cdot \1_{(0,\infty)}(h)\]
GeoStatsFunctions.CubicVariogram — TypeCubicVariogram(; range, sill, nugget)
+CubicVariogram(ball; sill, nugget)A cubic variogram with range, sill and nugget. Optionally, use a custom metric ball.
varioplot(CubicVariogram())
\[\gamma(h) = (s - n) \left[\left(\frac{15}{8}\left(\frac{h}{r}\right) - \frac{5}{4}\left(\frac{h}{r}\right)^3 + \frac{3}{8}\left(\frac{h}{r}\right)^5\right) \cdot \1_{(0,r)}(h) + \1_{[r,\infty)}(h)\right] + n \cdot \1_{(0,\infty)}(h)\]
GeoStatsFunctions.PentasphericalVariogram — TypePentasphericalVariogram(; range, sill, nugget)
+PentasphericalVariogram(ball; sill, nugget)A pentaspherical variogram with range, sill and nugget. Optionally, use a custom metric ball.
varioplot(PentasphericalVariogram())
\[\gamma(h) = (s - n) \left[1 - \frac{\sin(\pi h / r)}{\pi h / r}\right] + n \cdot \1_{(0,\infty)}(h)\]
GeoStatsFunctions.SineHoleVariogram — TypeSineHoleVariogram(; range, sill, nugget)
+SineHoleVariogram(ball; sill, nugget)A sine hole variogram with range, sill and nugget. Optionally, use a custom metric ball.
varioplot(SineHoleVariogram())
\[\gamma(h) = (s - n) \left[\left(1 - \frac{2}{\pi} \cos^{-1}\left(\frac{h}{r}\right) + \frac{2h}{\pi r} \sqrt{1 - \frac{h^2}{r^2}} \right) \cdot \1_{(0,r)}(h) + \1_{[r,\infty)}(h)\right] + n \cdot \1_{(0,\infty)}(h)\]
GeoStatsFunctions.CircularVariogram — TypeCircularVariogram(; range, sill, nugget)
+CircularVariogram(ball; sill, nugget)A circular variogram with range, sill and nugget. Optionally, use a custom metric ball.
varioplot(CircularVariogram())
\[\gamma(h) = sh^a + n \cdot \1_{(0,\infty)}(h)\]
GeoStatsFunctions.PowerVariogram — TypePowerVariogram(; scaling, exponent, nugget)A power variogram with scaling, exponent and nugget.
varioplot(PowerVariogram())
\[\gamma(h) = n \cdot \1_{(0,\infty)}(h)\]
GeoStatsFunctions.NuggetEffect — TypeNuggetEffect(; nugget)A pure nugget effect variogram with nugget.
varioplot(NuggetEffect(1.0))
Anisotropic models are easily obtained by defining an ellipsoid metric in place of the default Euclidean metric as shown in the following example. First, we create an ellipsoid that specifies the ranges and angles of rotation:
ellipsoid = MetricBall((3.0, 2.0, 1.0), RotZXZ(0.0, 0.0, 0.0))MetricBall((3.0 m, 2.0 m, 1.0 m), Mahalanobis)All rotations from Rotations.jl are supported as well as the following additional rotations from commercial geostatistical software:
GeoStatsBase.MinesightAngles — TypeMinesightAngles(θ₁, θ₂, θ₃)MineSight z-x'-y'' intrinsic rotation convention following the right-hand rule. All angles are in degrees and the sign convention is CW, CCW, CW positive.
The first rotation θ₁ is a horizontal rotation around the Z-axis, with positive being clockwise. The second rotation θ₂ is a rotation around the new X-axis, which moves the Y-axis into the desired position, with positive direction of rotation is up. The third rotation θ₃ is a rotation around the new Y-axis, which moves the X-axis into the desired position, with positive direction of rotation is up.
References
GeoStatsBase.DatamineAngles — TypeDatamineAngles(θ₁, θ₂, θ₃)Datamine ZXY rotation convention following the left-hand rule. All angles are in degrees and the signal convention is CW, CW, CW positive. Y is the principal axis.
GeoStatsBase.VulcanAngles — TypeVulcanAngles(θ₁, θ₂, θ₃)Vulcan ZYX rotation convention following the right-hand rule. All angles are in degrees and the signal convention is CW, CCW, CW positive. X is the principal axis.
GeoStatsBase.GslibAngles — TypeGslibAngles(θ₁, θ₂, θ₃)GSLIB z-x'-y'' intrinsic rotation convention following the left-hand rule. All angles are in degrees and the sign convention is CW, CCW, CW positive. Y is the principal axis.
The first rotation θ₁ is a rotation around the Z-axis, this is also called the azimuth. The second rotation θ₂ is a rotation around the new X-axis, this is also called the dip. The third rotation θ₃ is a rotation around the new Y-axis, this is also called the tilt.
References
We pass the ellipsoid as the first argument to the variogram model instead of specifying a single range with a keyword argument:
GaussianVariogram(ellipsoid, sill = 2.0)GaussianVariogram
+├─ sill: 2.0
+├─ nugget: 0.0
+├─ ranges: (3.0 m, 2.0 m, 1.0 m)
+└─ distance: MahalanobisTo illustrate the concept, consider the following 2D data set:
points = rand(Point, 50, crs=Cartesian2D) |> Scale(100)
+geotable = georef((Z=rand(50),), points)
+
+viz(geotable.geometry, color = geotable.Z)
We interpolate the data with different ellipsoids by varying the angle of rotation from $0$ to $2\pi$ clockwise:
θs = range(0.0, step = π/4, stop = 2π - π/4)
+
+# initialize figure
+fig = Mke.Figure(size = (800, 1600))
+
+# helper function to position subfigures
+pos = i -> CartesianIndices((4, 2))[i].I
+
+# Kriging with different angles
+for (i, θ) in enumerate(θs)
+ # ellipsoid rotated clockwise by angle θ
+ e = MetricBall((20.,5.), Angle2d(θ))
+
+ # anisotropic variogram model
+ γ = GaussianVariogram(e)
+
+ # domain of interpolation
+ grid = CartesianGrid(100, 100)
+
+ # perform interpolation
+ sol = geotable |> Interpolate(grid, Kriging(γ))
+
+ # visualize solution at subplot i
+ viz(fig[pos(i)...], sol.geometry, color = sol.Z,
+ axis = (title = "θ = $(rad2deg(θ))ᵒ",)
+ )
+end
+
+fig
A nested variogram model $\gamma(h) = c_1\cdot\gamma_1(h) + c_2\cdot\gamma_2(h) + \cdots + c_n\cdot\gamma_n(h)$ can be constructed from multiple variogram models, including matrix coefficients. The individual structures can be recovered in canonical form with the structures function:
GeoStatsFunctions.structures — Functionstructures(γ)Return the individual structures of a (possibly nested) variogram as a tuple. The structures are the total nugget cₒ, and the coefficients (or contributions) c[i] for the remaining non-trivial structures g[i] after normalization (i.e. sill=1, nugget=0).
Examples
γ₁ = GaussianVariogram(nugget=1, sill=2)
+γ₂ = SphericalVariogram(nugget=2, sill=3)
+
+γ = 2γ₁ + 3γ₂
+
+cₒ, c, g = structures(γ)γ₁ = GaussianVariogram(nugget=1, sill=2)
+γ₂ = ExponentialVariogram(nugget=2, sill=3)
+
+# nested model with matrix coefficients
+γ = [1.0 0.0; 0.0 1.0] * γ₁ + [2.0 0.5; 0.5 3.0] * γ₂
+
+# structures in canonical form
+cₒ, c, g = structures(γ)
+
+cₒ # matrix nugget2×2 Matrix{Float64}:
+ 5.0 1.0
+ 1.0 7.0c # matrix coefficients([1.0 0.0; 0.0 1.0], [2.0 0.5; 0.5 3.0])g # normalized structures(GaussianVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean), ExponentialVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean))Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.
The framework provides powerful visualization recipes for geospatial data science via the Makie.jl project. These recipes were carefully designed to maximize productivity and to protect users from GIS jargon. The main entry point is the viz function:
Meshes.viz — Functionviz(object; [options])Visualize Meshes.jl object with various options.
Available options
color - color of geometriesalpha - transparency in [0,1]colormap - color scheme/map from ColorSchemes.jlcolorrange - minimum and maximum color valuesshowsegments - visualize segmentssegmentcolor - color of segmentssegmentsize - width of segmentsshowpoints - visualize pointspointmarker - marker of pointspointcolor - color of pointspointsize - size of pointsThe option color can be a single scalar or a vector of scalars. For Mesh subtypes, the length of the vector of colors determines if the colors should be assigned to vertices or to elements.
Examples
Different coloring methods (vertex vs. element):
# vertex coloring (i.e. linear interpolation)
+viz(mesh, color = 1:nvertices(mesh))
+
+# element coloring (i.e. discrete colors)
+viz(mesh, color = 1:nelements(mesh))Different strategies to show the boundary of geometries (showsegments vs. boundary):
# visualize boundary with showsegments
+viz(polygon, showsegments = true)
+
+# visualize boundary with separate call
+viz(polygon)
+viz!(boundary(polygon))Notes
Meshes.viz! — Functionviz!(object; [options])Visualize Meshes.jl object in an existing scene with options forwarded to viz.
This function takes a geospatial domain as input and provides a set of aesthetic options to style the elements (i.e. geometries) of the domain.
Notice that the geometry column of our geospatial data type is a domain (i.e. data.geometry isa Domain), and that this design enables several optimizations in the visualization itself.
Users can also call Makie's plot function in the geometry column as in
Mke.plot(data.geometry)and this is equivalent to calling the viz recipe above. The plot function also works with various other objects such as EmpiricalHistogram and EmpiricalVariogram. That is convenient if you don't remember the name of the recipe.
Additionaly, we provide a basic scientific viewer to visualize all viewable variables in the data:
GeoTables.viewer — Functionviewer(geotable; kwargs...)Basic scientific viewer for geospatial table geotable.
Aesthetic options are forwarded via kwargs to the Meshes.viz recipe.
Other plots are listed below that can be useful for geostatistical analysis.
A hscatter plot between two variables var1 and var2 (possibly with var2 = var1) is a simple scatter plot in which the dots represent all ordered pairs of values of var1 and var2 at a given lag h.
𝒟 = 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
The PairPlots.jl package provides the pairplot function that can be used with any table, including tables of attributes obtained with the values function.
The Biplot.jl package provides 2D and 3D statistical biplots.
Settings
This document was generated with Documenter.jl version 1.8.0 on Thursday 14 November 2024. Using Julia version 1.11.1.