geossm.covmodel.spdeAppoxCov

class geossm.covmodel.spdeAppoxCov(latlon=True, geo_scale=57.29577951308232, nu=1, var=1.0, rescale=1.0, verbose=True)[source]

Bases: Matern

The SPDE approximation of the Matérn covariance model.

See [Rasmussen2003] for background on the Matérn family of covariance functions in the Gaussian-process setting that this class approximates. Solves the SPDE \((\kappa^2 - \Delta)^{\alpha/2} x = W\) by a finite element (FEM) discretisation with Neumann boundary conditions on a mesh covering the domain (see setup), giving the sparse precision matrix precision() – optionally marginalised onto the “inner” vertices of domain via precision(marginal=True), which is what other models (e.g. LRStateSpaceModel) should generally use.

Notes

spdeAppoxCov extends gstools.covmodel.Matern, but only precision() (and emp_range, sigma2k, distance()) are re-defined in terms of the FEM discretisation. The methods inherited unchanged from Matern – variogram, covariance, correlation, spectrum, spectral_density, and the *_yadrenko variants – describe the nominal, stationary, infinite-domain Matérn process that this SPDE approximates, using rescale as \(\kappa\). They are exact only in the far interior of the mesh, away from the boundary; they do not reflect the finite mesh’s boundary effects, nor the Schur-complement marginalisation used by precision(marginal=True). Use them for diagnostics/plotting against the nominal target kernel, not as the exact model – for that, use precision()/precision(marginal=True).

References

[Rasmussen2003]

Rasmussen, C. E., “Gaussian processes in machine learning.” Summer school on machine learning. Springer, Berlin, Heidelberg, (2003)

Parameters:
  • dim (int, optional) – dimension of the model. Includes the temporal dimension if temporal is true. To specify only the spatial dimension in that case, use spatial_dim. Default: 3

  • var (float, optional) – variance of the model (the nugget is not included in “this” variance) Default: 1.0

  • len_scale (float or list, optional) – length scale of the model. If a single value is given, the same length-scale will be used for every direction. If multiple values (for main and transversal directions) are given, anis will be recalculated accordingly. If only two values are given in 3D, the latter one will be used for both transversal directions. Default: 1.0

  • nugget (float, optional) – nugget of the model. Default: 0.0

  • anis (float or list, optional) –

    anisotropy ratios in the transversal directions [e_y, e_z].

    • e_y = l_y / l_x

    • e_z = l_z / l_x

    If only one value is given in 3D, e_y will be set to 1. This value will be ignored, if multiple len_scales are given. Default: 1.0

  • angles (float or list, optional) –

    angles of rotation (given in rad):

    • in 2D: given as rotation around z-axis

    • in 3D: given by yaw, pitch, and roll (known as Tait–Bryan angles)

    Default: 0.0

  • integral_scale (float or list or None, optional) – If given, len_scale will be ignored and recalculated, so that the integral scale of the model matches the given one. Default: None

  • rescale (float or None, optional) – Optional rescaling factor to divide the length scale with. This could be used for unit conversion or rescaling the length scale to coincide with e.g. the integral scale. Will be set by each model individually. Default: None

  • latlon (bool, optional) – Whether the model is describing 2D fields on earths surface described by latitude and longitude. When using this, the model will internally use the associated ‘Yadrenko’ model to represent a valid model. This means, the spatial distance \(r\) will be replaced by \(2\sin(\alpha/2)\), where \(\alpha\) is the great-circle distance, which is equal to the spatial distance of two points in 3D. As a consequence, dim will be set to 3 and anisotropy will be disabled. geo_scale can be set to e.g. earth’s radius, to have a meaningful len_scale parameter. Default: False

  • geo_scale (float, optional) – Geographic unit scaling in case of latlon coordinates to get a meaningful length scale unit. By default, len_scale is assumed to be in radians with latlon=True. Can be set to KM_SCALE to have len_scale in km or DEGREE_SCALE to have len_scale in degrees. Default: RADIAN_SCALE

  • temporal (bool, optional) – Create a metric spatio-temporal covariance model. Setting this to true will increase dim and field_dim by 1. spatial_dim will be field_dim - 1. The time-dimension is appended, meaning the pos tuple is (x,y,z,…,t). Default: False

  • spatial_dim (int, optional) – spatial dimension of the model. If given, the model dimension will be determined from this spatial dimension and the possible temporal dimension if temporal is ture. Default: None

  • hankel_kw (dict or None, optional) – Modify the init-arguments of hankel.SymmetricFourierTransform used for the spectrum calculation. Use with caution (Better: Don’t!). None is equivalent to {"a": -1, "b": 1, "N": 1000, "h": 0.001}. Default: None

  • **opt_arg – Optional arguments are covered by these keyword arguments. If present, they are described in the section Other Parameters.

Methods

covers_domain

Check that the mesh (set via setup()) covers domain_hull.

distance

Pairwise distance between mesh vertices (or points and mesh vertices).

generate_summary

Build the (label, value) rows shown by summary, as two lists of (str, [str]) pairs (left and right columns); includes the fem_solver's own rows once setup() has been called.

precision

Compute the precision matrix (Q) for the spatial process z with marginal variance sigma2_process.

print_info

Print msg prefixed with a UTC timestamp.

setup

Initialize the covariance model with a mesh.

summary

Return a statsmodels.iolib.summary.Summary table of the model (kappa/range/variance/nu) and, once setup() has been called, the underlying fem_solver's mesh diagnostics; also backs __str__.

Attributes

domain

Region(s) of scientific interest (see _validate_domain).

domain_hull

Convex hull of the union of domain -- a reasonable default choice of FEM computational domain Ω, e.g. to pass as domain to buildMesh2d before calling setup().

emp_range

Practical spatial correlation range of the Matérn field, \(\sqrt{8 \nu} / \kappa\), i.e. the same rescale-based range formula used by the parent Matern model, expressed in the distance units of rescale/geo_scale.

fem_solver

The FEMSolver built by setup from meshIO/domain.

meshIO

The meshio.Mesh passed to setup.

sigma2k

Marginal variance of the standardized SPDE approximation, i.e. the variance of the approximate field x(u) obtained from rescale (\(\kappa\)) alone.

property meshIO

The meshio.Mesh passed to setup.

Raises:

RuntimeError – If setup hasn’t been called yet.

property fem_solver

The FEMSolver built by setup from meshIO/domain.

Raises:

RuntimeError – If setup hasn’t been called yet.

setup(mesh_obj, domain=None, stats=True)[source]

Initialize the covariance model with a mesh.

This method must be called after instantiation to set up the FEM discretization.

Parameters:
  • mesh_obj (meshio.Mesh) – A meshio.Mesh object (already loaded in memory)

  • domain (shapely.geometry.Polygon or MultiPolygon, or a list/tuple) – of them, optional Region(s) of scientific interest, e.g. one polygon (or a single MultiPolygon) per landmass for a multi-polygon country (Sicily, Sardinia, mainland Italy). Any MultiPolygon is expanded into its constituent polygons. Forwarded to FEMSolver, where it is used only to classify mesh vertices as inner/outer – see _validate_domain for the full contract, including how it differs from the FEM computational domain \(\Omega\) (implicitly given by mesh_obj; see domain_hull for a reasonable default choice used to build it). If omitted, FEMSolver falls back to the convex hull of the mesh’s own vertices, i.e. every vertex is treated as inner.

  • stats (bool, optional) – Whether to compute mesh quality statistics (default is True)

Returns:

self – Returns self for method chaining

Raises:
  • ValueError – If domain is invalid (see _validate_domain)

  • IOError – If meshpath file cannot be read

  • RuntimeError – If mesh is invalid or FEM setup fails

Examples

>>> mesh = meshio.read("my_mesh.msh")
>>> cov.setup(mesh_obj=mesh)
precision(rescale=None, marginal=False)[source]

Compute the precision matrix (Q) for the spatial process z with marginal variance sigma2_process.

Parameters:
  • rescale (float, optional) – If given, updates self.rescale (the SPDE kappa) before computing the precision matrix.

  • marginal (bool, optional) – If True, return the precision matrix marginalised (via the Schur complement, see _schur_marginal_precision) onto the domain inner vertices only, integrating out the outer (boundary) vertices – i.e. the actual precision of the latent field once the mesh’s boundary effects have been removed, as used e.g. by LRStateSpaceModel. Default: False, i.e. the full FEM precision matrix over every mesh vertex (inner and outer).

Returns:

scipy.sparse matrix if `marginal=False`, else a dense numpy array – (the Schur complement is generally dense even though Q is sparse).

distance(points=None)[source]

Pairwise distance between mesh vertices (or points and mesh vertices).

Honors latlon/geo_scale the same way the parent Matern model does: great-circle distance (in geo_scale units) when latlon=True, flat Euclidean distance on the raw mesh/point coordinates otherwise.

Note this only affects distance-based diagnostics (e.g. comparing against variogram/covariance, see the class Notes) – the FEM stiffness/mass matrices behind precision() are always assembled directly on the raw mesh coordinates with a flat metric, regardless of latlon. A raw lon/lat mesh under latlon=True therefore still carries anisotropic distortion (up to ~cos(latitude)) in the FEM operator that this distance correction does not fix; for a geodesically accurate SPDE solve, build the mesh in a projected (e.g. UTM/km) CRS instead.

property emp_range

Practical spatial correlation range of the Matérn field, \(\sqrt{8 \nu} / \kappa\), i.e. the same rescale-based range formula used by the parent Matern model, expressed in the distance units of rescale/geo_scale.

property sigma2k

Marginal variance of the standardized SPDE approximation, i.e. the variance of the approximate field x(u) obtained from rescale (\(\kappa\)) alone. Used by _compute_precision_spde as the scalar multiplier of the precision matrix.

property domain

Region(s) of scientific interest (see _validate_domain). Delegates to fem_solver.domain, so it is only available once setup() has been called – domain is now a setup()-only argument, since nothing reads it before the mesh exists.

property domain_hull

Convex hull of the union of domain – a reasonable default choice of FEM computational domain Ω, e.g. to pass as domain to buildMesh2d before calling setup(). See _domain_hull.

Raises:

RuntimeError – If setup() hasn’t been called yet – domain is unknown until then.

covers_domain()[source]

Check that the mesh (set via setup()) covers domain_hull. See FEMSolver.covers_domain.

generate_summary()[source]

Build the (label, value) rows shown by summary, as two lists of (str, [str]) pairs (left and right columns); includes the fem_solver’s own rows once setup() has been called.

summary()[source]

Return a statsmodels.iolib.summary.Summary table of the model (kappa/range/variance/nu) and, once setup() has been called, the underlying fem_solver’s mesh diagnostics; also backs __str__.

print_info(msg)[source]

Print msg prefixed with a UTC timestamp.