16463739022

Committed 23 Jul 2025 06:57AM UTC coverage: 93.652%. Remained the same

Build # 16463739022

Build Type

Pull #389

github

Committed by

web-flow

Commit Message

Merge 956658981 into 6d179b230

Pull Request Pull Request #389: Fully support Python version 3.8 again

Coverage Stats

12 of 24 new or added lines in 6 files covered. (50.0%)

3659 of 3907 relevant lines covered (93.65%)

0.94 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

97.26

/src/gstools/variogram/variogram.py

"""
GStools subpackage providing tools for estimating and fitting variograms.

.. currentmodule:: gstools.variogram.variogram

The following functions are provided

.. autosummary::
   vario_estimate
   vario_estimate_axis
"""

import numpy as np
from gstools_cython.variogram import directional as directional_c
from gstools_cython.variogram import ma_structured as ma_structured_c
from gstools_cython.variogram import structured as structured_c
from gstools_cython.variogram import unstructured as unstructured_c

from gstools import config
from gstools.normalizer.tools import remove_trend_norm_mean
from gstools.tools import RADIAN_SCALE
from gstools.tools.geometric import (
    ang2dir,
    format_struct_pos_shape,
    format_unstruct_pos_shape,
    generate_grid,
)
from gstools.variogram.binning import standard_bins

if config._GSTOOLS_CORE_AVAIL:  # pragma: no cover
    from gstools_core import variogram_directional as directional_gsc
    from gstools_core import variogram_ma_structured as ma_structured_gsc
    from gstools_core import variogram_structured as structured_gsc
    from gstools_core import variogram_unstructured as unstructured_gsc

__all__ = [
    "vario_estimate",
    "vario_estimate_axis",
    "vario_estimate_unstructured",
    "vario_estimate_structured",
]


AXIS = ["x", "y", "z"]
AXIS_DIR = {"x": 0, "y": 1, "z": 2}


def _directional(
    field,
    bin_edges,
    pos,
    direction,
    angles_tol=np.pi / 8.0,
    bandwidth=-1.0,
    separate_dirs=False,
    estimator_type="m",
    num_threads=None,
):
    """A wrapper function for calling the directional variogram algorithms."""
    if config.USE_GSTOOLS_CORE and config._GSTOOLS_CORE_AVAIL:
        directional_fct = directional_gsc
    else:
        directional_fct = directional_c
    return directional_fct(
        field,
        bin_edges,
        pos,
        direction,
        angles_tol,
        bandwidth,
        separate_dirs,
        estimator_type,
        num_threads,
    )


def _unstructured(
    field,
    bin_edges,
    pos,
    estimator_type="m",
    distance_type="e",
    num_threads=None,
):
    """A wrapper function for calling the unstructured variogram algorithms."""
    if config.USE_GSTOOLS_CORE and config._GSTOOLS_CORE_AVAIL:
        unstructured_fct = unstructured_gsc
    else:
        unstructured_fct = unstructured_c
    return unstructured_fct(
        field,
        bin_edges,
        pos,
        estimator_type,
        distance_type,
        num_threads,
    )


def _structured(
    field,
    estimator_type="m",
    num_threads=None,
):
    """A wrapper function for calling the structured variogram algorithms."""
    if config.USE_GSTOOLS_CORE and config._GSTOOLS_CORE_AVAIL:
        structured_fct = structured_gsc
    else:
        structured_fct = structured_c
    return structured_fct(field, estimator_type, num_threads)


def _ma_structured(
    field,
    mask,
    estimator_type="m",
    num_threads=None,
):
    """A wrapper function for calling the masked struct. variogram algorithms."""
    if config.USE_GSTOOLS_CORE and config._GSTOOLS_CORE_AVAIL:
        ma_structured_fct = ma_structured_gsc
    else:
        ma_structured_fct = ma_structured_c
    return ma_structured_fct(field, mask, estimator_type, num_threads)


def _set_estimator(estimator):
    """Translate the verbose Python estimator identifier to single char."""
    if estimator.lower() == "matheron":
        cython_estimator = "m"
    elif estimator.lower() == "cressie":
        cython_estimator = "c"
    else:
        raise ValueError(f"Unknown variogram estimator function: {estimator}")
    return cython_estimator


def _separate_dirs_test(direction, angles_tol):
    """Check if given directions are separated."""
    if direction is None or direction.shape[0] < 2:
        return True
    separate_dirs = True
    for i in range(direction.shape[0] - 1):
        for j in range(i + 1, direction.shape[0]):
            s_prod = np.minimum(np.abs(np.dot(direction[i], direction[j])), 1)
            separate_dirs &= np.arccos(s_prod) >= 2 * angles_tol
    # gstools-core doesn't like the type `numpy.bool_`
    return bool(separate_dirs)


def vario_estimate(
    pos,
    field,
    bin_edges=None,
    sampling_size=None,
    sampling_seed=None,
    estimator="matheron",
    latlon=False,
    direction=None,
    angles=None,
    angles_tol=np.pi / 8,
    bandwidth=None,
    no_data=np.nan,
    mask=np.ma.nomask,
    mesh_type="unstructured",
    return_counts=False,
    mean=None,
    normalizer=None,
    trend=None,
    fit_normalizer=False,
    geo_scale=RADIAN_SCALE,
    **std_bins,
):
    r"""
    Estimates the empirical variogram.

    The algorithm calculates following equation:

    .. math::
       \gamma(r_k) = \frac{1}{2 N(r_k)} \sum_{i=1}^{N(r_k)} (z(\mathbf x_i) -
       z(\mathbf x_i'))^2 \; ,

    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
    being the bins.

    Or if the estimator "cressie" was chosen:

    .. math::
       \gamma(r_k) = \frac{\frac{1}{2}\left(\frac{1}{N(r_k)}\sum_{i=1}^{N(r_k)}
       \left|z(\mathbf x_i) - z(\mathbf x_i')\right|^{0.5}\right)^4}
       {0.457 + 0.494 / N(r_k) + 0.045 / N^2(r_k)} \; ,

    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
    being the bins.
    The Cressie estimator is more robust to outliers [Webster2007]_.

    By providing `direction` vector[s] or angles, a directional variogram
    can be calculated. If multiple directions are given, a set of variograms
    will be returned.
    Directional bining is controlled by a given angle tolerance (`angles_tol`)
    and an optional `bandwidth`, that truncates the width of the search band
    around the given direction[s].

    To reduce the calculation time, `sampling_size` could be passed to sample
    down the number of field points.

    Parameters
    ----------
    pos : :class:`list`
        the position tuple, containing either the point coordinates (x, y, ...)
        or the axes descriptions (for mesh_type='structured')
    field : :class:`numpy.ndarray` or :class:`list` of :class:`numpy.ndarray`
        The spatially distributed data.
        Can also be of type :class:`numpy.ma.MaskedArray` to use masked values.
        You can pass a list of fields, that will be used simultaneously.
        This could be helpful, when there are multiple realizations at the
        same points, with the same statistical properties.
    bin_edges : :class:`numpy.ndarray`, optional
        the bins on which the variogram will be calculated.
        If :any:`None` are given, standard bins provided by the
        :any:`standard_bins` routine will be used. Default: :any:`None`
    sampling_size : :class:`int` or :any:`None`, optional
        for large input data, this method can take a long
        time to compute the variogram, therefore this argument specifies
        the number of data points to sample randomly
        Default: :any:`None`
    sampling_seed : :class:`int` or :any:`None`, optional
        seed for samples if sampling_size is given.
        Default: :any:`None`
    estimator : :class:`str`, optional
        the estimator function, possible choices:

            * "matheron": the standard method of moments of Matheron
            * "cressie": an estimator more robust to outliers

        Default: "matheron"
    latlon : :class:`bool`, optional
        Whether the data is representing 2D fields on earths surface described
        by latitude and longitude. When using this, the estimator will
        use great-circle distance for variogram estimation.
        Note, that only an isotropic variogram can be estimated and a
        ValueError will be raised, if a direction was specified.
        Bin edges need to be given in radians in this case.
        Default: False
    direction : :class:`list` of :class:`numpy.ndarray`, optional
        directions to evaluate a directional variogram.
        Angular tolerance is given by `angles_tol`.
        bandwidth to cut off how wide the search for point pairs should be
        is given by `bandwidth`.
        You can provide multiple directions at once to get one variogram
        for each direction.
        For a single direction you can also use the `angles` parameter,
        to provide the direction by its spherical coordinates.
        Default: :any:`None`
    angles : :class:`numpy.ndarray`, optional
        the angles of the main axis to calculate the variogram for in radians
        angle definitions from ISO standard 80000-2:2009
        for 1d this parameter will have no effect at all
        for 2d supply one angle which is
        azimuth :math:`\varphi` (ccw from +x in xy plane)
        for 3d supply two angles which are
        azimuth :math:`\varphi` (ccw from +x in xy plane)
        and inclination :math:`\theta` (cw from +z).
        Can be used instead of direction.
        Default: :any:`None`
    angles_tol : class:`float`, optional
        the tolerance around the variogram angle to count a point as being
        within this direction from another point (the angular tolerance around
        the directional vector given by angles)
        Default: `np.pi/8` = 22.5°
    bandwidth : class:`float`, optional
        bandwidth to cut off the angular tolerance for directional variograms.
        If None is given, only the `angles_tol` parameter will control the
        point selection.
        Default: :any:`None`
    no_data : :class:`float`, optional
        Value to identify missing data in the given field.
        Default: `numpy.nan`
    mask : :class:`numpy.ndarray` of :class:`bool`, optional
        Mask to deselect data in the given field.
        Default: :any:`numpy.ma.nomask`
    mesh_type : :class:`str`, optional
        'structured' / 'unstructured', indicates whether the pos tuple
        describes the axis or the point coordinates.
        Default: `'unstructured'`
    return_counts: :class:`bool`, optional
        if set to true, this function will also return the number of data
        points found at each lag distance as a third return value
        Default: False
    mean : :class:`float`, optional
        mean value used to shift normalized input data.
        Can also be a callable. The default is None.
    normalizer : :any:`None` or :any:`Normalizer`, optional
        Normalizer to be applied to the input data to gain normality.
        The default is None.
    trend : :any:`None` or :class:`float` or :any:`callable`, optional
        A callable trend function. Should have the signature: f(x, [y, z, ...])
        If no normalizer is applied, this behaves equal to 'mean'.
        The default is None.
    fit_normalizer : :class:`bool`, optional
        Whether to fit the data-normalizer to the given (detrended) field.
        Default: False
    geo_scale : :class:`float`, optional
        Geographic unit scaling in case of latlon coordinates to get a
        meaningful bins unit.
        By default, bins are assumed to be in radians with latlon=True.
        Can be set to :any:`KM_SCALE` to have bins in km or
        :any:`DEGREE_SCALE` to have bins in degrees.
        Default: :any:`RADIAN_SCALE`
    **std_bins
        Optional arguments that are forwarded to the :any:`standard_bins` routine
        if no bins are given (bin_no, max_dist).

    Returns
    -------
    bin_centers : (n), :class:`numpy.ndarray`
        The bin centers.
    gamma : (n) or (d, n), :class:`numpy.ndarray`
        The estimated variogram values at bin centers.
        Is stacked if multiple `directions` (d>1) are given.
    counts : (n) or (d, n), :class:`numpy.ndarray`, optional
        The number of point pairs found for each bin.
        Is stacked if multiple `directions` (d>1) are given.
        Only provided if `return_counts` is True.
    normalizer : :any:`Normalizer`, optional
        The fitted normalizer for the given data.
        Only provided if `fit_normalizer` is True.

    Notes
    -----
    Internally uses double precision and also returns doubles.

    References
    ----------
    .. [Webster2007] Webster, R. and Oliver, M. A.
           "Geostatistics for environmental scientists.",
           John Wiley & Sons. (2007)
    """
    if bin_edges is not None:
        bin_edges = np.atleast_1d(np.asarray(bin_edges, dtype=np.double))
        bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0
    # allow multiple fields at same positions (ndmin=2: first axis -> field ID)
    # need to convert to ma.array, since list of ma.array is not recognised
    field = np.ma.array(field, ndmin=2, dtype=np.double, copy=True)
    masked = np.ma.is_masked(field) or np.any(mask)
    # catch special case if everything is masked
    if masked and np.all(mask):
        bin_centers = np.empty(0) if bin_edges is None else bin_centers
        estimates = np.zeros_like(bin_centers)
        if return_counts:
            return bin_centers, estimates, np.zeros_like(estimates, dtype=int)
        return bin_centers, estimates
    if not masked:
        field = field.filled()
    # check mesh shape
    if mesh_type != "unstructured":
        pos, __, dim = format_struct_pos_shape(
            pos, field.shape, check_stacked_shape=True
        )
        pos = generate_grid(pos)
    else:
        pos, __, dim = format_unstruct_pos_shape(
            pos, field.shape, check_stacked_shape=True
        )
    if latlon and dim != 2:
        raise ValueError("Variogram: given field needs to be 2D for lat-lon.")
    # prepare the field
    pnt_cnt = len(pos[0])
    field = field.reshape((-1, pnt_cnt))
    # apply mask if wanted
    if masked:
        # if fields have different masks, take the minimal common mask
        # given mask will be applied in addition
        # selected region is the inverted masked (unmasked values)
        if np.size(mask) > 1:  # not only np.ma.nomask
            select = np.invert(
                np.logical_or(
                    np.reshape(mask, pnt_cnt), np.all(field.mask, axis=0)
                )
            )
        else:
            select = np.invert(np.all(field.mask, axis=0))
        pos = pos[:, select]
        field.fill_value = np.nan  # use no-data val. for remaining masked vals
        field = field[:, select].filled()  # convert to ndarray
        select = mask = None  # free space
        pnt_cnt = len(pos[0])  # pnt cnt reduced now
    # set no_data values
    if not np.isnan(no_data):
        field[np.isclose(field, float(no_data))] = np.nan
    # set directions
    dir_no = 0
    if direction is not None and dim > 1:
        direction = np.atleast_2d(np.asarray(direction, dtype=np.double))
        if len(direction.shape) > 2:
            raise ValueError(f"Can't interpret directions: {direction}")
        if direction.shape[1] != dim:
            raise ValueError(f"Can't interpret directions: {direction}")
        dir_no = direction.shape[0]
    # convert given angles to direction vector
    if angles is not None and direction is None and dim > 1:
        direction = ang2dir(angles=angles, dtype=np.double, dim=dim)
        dir_no = direction.shape[0]
    # prepare directional variogram
    if dir_no > 0:
        if latlon:
            raise ValueError("Directional variogram not allowed for lat-lon.")
        norms = np.linalg.norm(direction, axis=1)
        if np.any(np.isclose(norms, 0)):
            raise ValueError(f"Zero length directions: {direction}")
        # only unit-vectors for directions
        direction = np.divide(direction, norms[:, np.newaxis])
        # negative bandwidth to turn it off
        bandwidth = float(bandwidth) if bandwidth is not None else -1.0
        angles_tol = float(angles_tol)
    # prepare sampled variogram
    if sampling_size is not None and sampling_size < pnt_cnt:
        sampled_idx = np.random.RandomState(sampling_seed).choice(
            np.arange(pnt_cnt), sampling_size, replace=False
        )
        field = field[:, sampled_idx]
        pos = pos[:, sampled_idx]
    # create bins
    if bin_edges is None:
        bin_edges = standard_bins(
            pos, dim, latlon, geo_scale=geo_scale, **std_bins
        )
        bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0
    if latlon:
        # internally we always use radians
        bin_edges /= geo_scale
    # normalize field
    norm_field_out = remove_trend_norm_mean(
        *(pos, field, mean, normalizer, trend),
        check_shape=False,
        stacked=True,
        fit_normalizer=fit_normalizer,
    )
    field = norm_field_out[0] if fit_normalizer else norm_field_out
    norm_out = (norm_field_out[1],) if fit_normalizer else ()
    # select variogram estimator
    cython_estimator = _set_estimator(estimator)
    # run
    if dir_no == 0:
        # "h"aversine or "e"uclidean distance type
        distance_type = "h" if latlon else "e"
        estimates, counts = _unstructured(
            field,
            bin_edges,
            pos,
            estimator_type=cython_estimator,
            distance_type=distance_type,
            num_threads=config.NUM_THREADS,
        )
    else:
        estimates, counts = _directional(
            field,
            bin_edges,
            pos,
            direction,
            angles_tol,
            bandwidth,
            separate_dirs=_separate_dirs_test(direction, angles_tol),
            estimator_type=cython_estimator,
            num_threads=config.NUM_THREADS,
        )
        if dir_no == 1:
            estimates, counts = estimates[0], counts[0]
    est_out = (estimates, counts)
    return (bin_centers,) + est_out[: 2 if return_counts else 1] + norm_out


def vario_estimate_axis(
    field, direction="x", estimator="matheron", no_data=np.nan
):
    r"""Estimates the variogram along array axis.

    The indices of the given direction are used for the bins.
    Uniform spacings along the given axis are assumed.

    The algorithm calculates following equation:

    .. math::
       \gamma(r_k) = \frac{1}{2 N(r_k)} \sum_{i=1}^{N(r_k)} (z(\mathbf x_i) -
       z(\mathbf x_i'))^2 \; ,

    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
    being the bins.

    Or if the estimator "cressie" was chosen:

    .. math::
       \gamma(r_k) = \frac{\frac{1}{2}\left(\frac{1}{N(r_k)}\sum_{i=1}^{N(r_k)}
       \left|z(\mathbf x_i) - z(\mathbf x_i')\right|^{0.5}\right)^4}
       {0.457 + 0.494 / N(r_k) + 0.045 / N^2(r_k)} \; ,

    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
    being the bins.
    The Cressie estimator is more robust to outliers [Webster2007]_.

    Parameters
    ----------
    field : :class:`numpy.ndarray` or :class:`numpy.ma.MaskedArray`
        the spatially distributed data (can be masked)
    direction : :class:`str` or :class:`int`
        the axis over which the variogram will be estimated (x, y, z)
        or (0, 1, 2, ...)
    estimator : :class:`str`, optional
        the estimator function, possible choices:

            * "matheron": the standard method of moments of Matheron
            * "cressie": an estimator more robust to outliers

        Default: "matheron"

    no_data : :class:`float`, optional
        Value to identify missing data in the given field.
        Default: `numpy.nan`

    Returns
    -------
    :class:`numpy.ndarray`
        the estimated variogram along the given direction.

    Warnings
    --------
    It is assumed that the field is defined on an equidistant Cartesian grid.

    Notes
    -----
    Internally uses double precision and also returns doubles.

    References
    ----------
    .. [Webster2007] Webster, R. and Oliver, M. A.
           "Geostatistics for environmental scientists.",
           John Wiley & Sons. (2007)
    """
    missing_mask = (
        np.isnan(field) if np.isnan(no_data) else np.isclose(field, no_data)
    )
    missing = np.any(missing_mask)
    masked = np.ma.is_masked(field) or missing
    if masked:
        field = np.ma.array(field, ndmin=1, dtype=np.double)
        if missing:
            field.mask = np.logical_or(field.mask, missing_mask)
        mask = np.ma.getmaskarray(field)
    else:
        field = np.atleast_1d(np.asarray(field, dtype=np.double))
        missing_mask = None  # free space

    axis_to_swap = AXIS_DIR[direction] if direction in AXIS else int(direction)
    # desired axis first, convert to 2D array afterwards
    field = field.swapaxes(0, axis_to_swap)
    field = field.reshape((field.shape[0], -1))
    if masked:
        mask = mask.swapaxes(0, axis_to_swap)
        mask = mask.reshape((mask.shape[0], -1))

    cython_estimator = _set_estimator(estimator)

    if masked:
        return _ma_structured(
            field, mask, cython_estimator, num_threads=config.NUM_THREADS
        )
    return _structured(field, cython_estimator, num_threads=config.NUM_THREADS)


# for backward compatibility
vario_estimate_unstructured = vario_estimate
vario_estimate_structured = vario_estimate_axis

1	"""
2	GStools subpackage providing tools for estimating and fitting variograms.
3
4	.. currentmodule:: gstools.variogram.variogram
5
6	The following functions are provided
7
8	.. autosummary::
9	vario_estimate
10	vario_estimate_axis
11	"""
12
13	import numpy as np	1✔
14	from gstools_cython.variogram import directional as directional_c	1✔
15	from gstools_cython.variogram import ma_structured as ma_structured_c	1✔
16	from gstools_cython.variogram import structured as structured_c	1✔
17	from gstools_cython.variogram import unstructured as unstructured_c	1✔
18
19	from gstools import config	1✔
20	from gstools.normalizer.tools import remove_trend_norm_mean	1✔
21	from gstools.tools import RADIAN_SCALE	1✔
22	from gstools.tools.geometric import (	1✔
23	ang2dir,
24	format_struct_pos_shape,
25	format_unstruct_pos_shape,
26	generate_grid,
27	)
28	from gstools.variogram.binning import standard_bins	1✔
29
30	if config._GSTOOLS_CORE_AVAIL: # pragma: no cover
31	from gstools_core import variogram_directional as directional_gsc
32	from gstools_core import variogram_ma_structured as ma_structured_gsc
33	from gstools_core import variogram_structured as structured_gsc
34	from gstools_core import variogram_unstructured as unstructured_gsc
35
36	__all__ = [	1✔
37	"vario_estimate",
38	"vario_estimate_axis",
39	"vario_estimate_unstructured",
40	"vario_estimate_structured",
41	]
42
43
44	AXIS = ["x", "y", "z"]	1✔
45	AXIS_DIR = {"x": 0, "y": 1, "z": 2}	1✔
46
47
48	def _directional(	1✔
49	field,
50	bin_edges,
51	pos,
52	direction,
53	angles_tol=np.pi / 8.0,
54	bandwidth=-1.0,
55	separate_dirs=False,
56	estimator_type="m",
57	num_threads=None,
58	):
59	"""A wrapper function for calling the directional variogram algorithms."""
60	if config.USE_GSTOOLS_CORE and config._GSTOOLS_CORE_AVAIL:	1✔
NEW 61	directional_fct = directional_gsc	×
62	else:
63	directional_fct = directional_c	1✔
64	return directional_fct(	1✔
65	field,
66	bin_edges,
67	pos,
68	direction,
69	angles_tol,
70	bandwidth,
71	separate_dirs,
72	estimator_type,
73	num_threads,
74	)
75
76
77	def _unstructured(	1✔
78	field,
79	bin_edges,
80	pos,
81	estimator_type="m",
82	distance_type="e",
83	num_threads=None,
84	):
85	"""A wrapper function for calling the unstructured variogram algorithms."""
86	if config.USE_GSTOOLS_CORE and config._GSTOOLS_CORE_AVAIL:	1✔
NEW 87	unstructured_fct = unstructured_gsc	×
88	else:
89	unstructured_fct = unstructured_c	1✔
90	return unstructured_fct(	1✔
91	field,
92	bin_edges,
93	pos,
94	estimator_type,
95	distance_type,
96	num_threads,
97	)
98
99
100	def _structured(	1✔
101	field,
102	estimator_type="m",
103	num_threads=None,
104	):
105	"""A wrapper function for calling the structured variogram algorithms."""
106	if config.USE_GSTOOLS_CORE and config._GSTOOLS_CORE_AVAIL:	1✔
NEW 107	structured_fct = structured_gsc	×
108	else:
109	structured_fct = structured_c	1✔
110	return structured_fct(field, estimator_type, num_threads)	1✔
111
112
113	def _ma_structured(	1✔
114	field,
115	mask,
116	estimator_type="m",
117	num_threads=None,
118	):
119	"""A wrapper function for calling the masked struct. variogram algorithms."""
120	if config.USE_GSTOOLS_CORE and config._GSTOOLS_CORE_AVAIL:	1✔
NEW 121	ma_structured_fct = ma_structured_gsc	×
122	else:
123	ma_structured_fct = ma_structured_c	1✔
124	return ma_structured_fct(field, mask, estimator_type, num_threads)	1✔
125
126
127	def _set_estimator(estimator):	1✔
128	"""Translate the verbose Python estimator identifier to single char."""
129	if estimator.lower() == "matheron":	1✔
130	cython_estimator = "m"	1✔
131	elif estimator.lower() == "cressie":	1✔
132	cython_estimator = "c"	1✔
133	else:
134	raise ValueError(f"Unknown variogram estimator function: {estimator}")	1✔
135	return cython_estimator	1✔
136
137
138	def _separate_dirs_test(direction, angles_tol):	1✔
139	"""Check if given directions are separated."""
140	if direction is None or direction.shape[0] < 2:	1✔
141	return True	1✔
142	separate_dirs = True	1✔
143	for i in range(direction.shape[0] - 1):	1✔
144	for j in range(i + 1, direction.shape[0]):	1✔
145	s_prod = np.minimum(np.abs(np.dot(direction[i], direction[j])), 1)	1✔
146	separate_dirs &= np.arccos(s_prod) >= 2 * angles_tol	1✔
147	# gstools-core doesn't like the type `numpy.bool_`
148	return bool(separate_dirs)	1✔
149
150
151	def vario_estimate(	1✔
152	pos,
153	field,
154	bin_edges=None,
155	sampling_size=None,
156	sampling_seed=None,
157	estimator="matheron",
158	latlon=False,
159	direction=None,
160	angles=None,
161	angles_tol=np.pi / 8,
162	bandwidth=None,
163	no_data=np.nan,
164	mask=np.ma.nomask,
165	mesh_type="unstructured",
166	return_counts=False,
167	mean=None,
168	normalizer=None,
169	trend=None,
170	fit_normalizer=False,
171	geo_scale=RADIAN_SCALE,
172	**std_bins,
173	):
174	r"""
175	Estimates the empirical variogram.
176
177	The algorithm calculates following equation:
178
179	.. math::
180	\gamma(r_k) = \frac{1}{2 N(r_k)} \sum_{i=1}^{N(r_k)} (z(\mathbf x_i) -
181	z(\mathbf x_i'))^2 \; ,
182
183	with :math:`r_k \leq \\| \mathbf x_i - \mathbf x_i' \\| < r_{k+1}`
184	being the bins.
185
186	Or if the estimator "cressie" was chosen:
187
188	.. math::
189	\gamma(r_k) = \frac{\frac{1}{2}\left(\frac{1}{N(r_k)}\sum_{i=1}^{N(r_k)}
190	\left\|z(\mathbf x_i) - z(\mathbf x_i')\right\|^{0.5}\right)^4}
191	{0.457 + 0.494 / N(r_k) + 0.045 / N^2(r_k)} \; ,
192
193	with :math:`r_k \leq \\| \mathbf x_i - \mathbf x_i' \\| < r_{k+1}`
194	being the bins.
195	The Cressie estimator is more robust to outliers [Webster2007]_.
196
197	By providing `direction` vector[s] or angles, a directional variogram
198	can be calculated. If multiple directions are given, a set of variograms
199	will be returned.
200	Directional bining is controlled by a given angle tolerance (`angles_tol`)
201	and an optional `bandwidth`, that truncates the width of the search band
202	around the given direction[s].
203
204	To reduce the calculation time, `sampling_size` could be passed to sample
205	down the number of field points.
206
207	Parameters
208	----------
209	pos : :class:`list`
210	the position tuple, containing either the point coordinates (x, y, ...)
211	or the axes descriptions (for mesh_type='structured')
212	field : :class:`numpy.ndarray` or :class:`list` of :class:`numpy.ndarray`
213	The spatially distributed data.
214	Can also be of type :class:`numpy.ma.MaskedArray` to use masked values.
215	You can pass a list of fields, that will be used simultaneously.
216	This could be helpful, when there are multiple realizations at the
217	same points, with the same statistical properties.
218	bin_edges : :class:`numpy.ndarray`, optional
219	the bins on which the variogram will be calculated.
220	If :any:`None` are given, standard bins provided by the
221	:any:`standard_bins` routine will be used. Default: :any:`None`
222	sampling_size : :class:`int` or :any:`None`, optional
223	for large input data, this method can take a long
224	time to compute the variogram, therefore this argument specifies
225	the number of data points to sample randomly
226	Default: :any:`None`
227	sampling_seed : :class:`int` or :any:`None`, optional
228	seed for samples if sampling_size is given.
229	Default: :any:`None`
230	estimator : :class:`str`, optional
231	the estimator function, possible choices:
232
233	* "matheron": the standard method of moments of Matheron
234	* "cressie": an estimator more robust to outliers
235
236	Default: "matheron"
237	latlon : :class:`bool`, optional
238	Whether the data is representing 2D fields on earths surface described
239	by latitude and longitude. When using this, the estimator will
240	use great-circle distance for variogram estimation.
241	Note, that only an isotropic variogram can be estimated and a
242	ValueError will be raised, if a direction was specified.
243	Bin edges need to be given in radians in this case.
244	Default: False
245	direction : :class:`list` of :class:`numpy.ndarray`, optional
246	directions to evaluate a directional variogram.
247	Angular tolerance is given by `angles_tol`.
248	bandwidth to cut off how wide the search for point pairs should be
249	is given by `bandwidth`.
250	You can provide multiple directions at once to get one variogram
251	for each direction.
252	For a single direction you can also use the `angles` parameter,
253	to provide the direction by its spherical coordinates.
254	Default: :any:`None`
255	angles : :class:`numpy.ndarray`, optional
256	the angles of the main axis to calculate the variogram for in radians
257	angle definitions from ISO standard 80000-2:2009
258	for 1d this parameter will have no effect at all
259	for 2d supply one angle which is
260	azimuth :math:`\varphi` (ccw from +x in xy plane)
261	for 3d supply two angles which are
262	azimuth :math:`\varphi` (ccw from +x in xy plane)
263	and inclination :math:`\theta` (cw from +z).
264	Can be used instead of direction.
265	Default: :any:`None`
266	angles_tol : class:`float`, optional
267	the tolerance around the variogram angle to count a point as being
268	within this direction from another point (the angular tolerance around
269	the directional vector given by angles)
270	Default: `np.pi/8` = 22.5°
271	bandwidth : class:`float`, optional
272	bandwidth to cut off the angular tolerance for directional variograms.
273	If None is given, only the `angles_tol` parameter will control the
274	point selection.
275	Default: :any:`None`
276	no_data : :class:`float`, optional
277	Value to identify missing data in the given field.
278	Default: `numpy.nan`
279	mask : :class:`numpy.ndarray` of :class:`bool`, optional
280	Mask to deselect data in the given field.
281	Default: :any:`numpy.ma.nomask`
282	mesh_type : :class:`str`, optional
283	'structured' / 'unstructured', indicates whether the pos tuple
284	describes the axis or the point coordinates.
285	Default: `'unstructured'`
286	return_counts: :class:`bool`, optional
287	if set to true, this function will also return the number of data
288	points found at each lag distance as a third return value
289	Default: False
290	mean : :class:`float`, optional
291	mean value used to shift normalized input data.
292	Can also be a callable. The default is None.
293	normalizer : :any:`None` or :any:`Normalizer`, optional
294	Normalizer to be applied to the input data to gain normality.
295	The default is None.
296	trend : :any:`None` or :class:`float` or :any:`callable`, optional
297	A callable trend function. Should have the signature: f(x, [y, z, ...])
298	If no normalizer is applied, this behaves equal to 'mean'.
299	The default is None.
300	fit_normalizer : :class:`bool`, optional
301	Whether to fit the data-normalizer to the given (detrended) field.
302	Default: False
303	geo_scale : :class:`float`, optional
304	Geographic unit scaling in case of latlon coordinates to get a
305	meaningful bins unit.
306	By default, bins are assumed to be in radians with latlon=True.
307	Can be set to :any:`KM_SCALE` to have bins in km or
308	:any:`DEGREE_SCALE` to have bins in degrees.
309	Default: :any:`RADIAN_SCALE`
310	**std_bins
311	Optional arguments that are forwarded to the :any:`standard_bins` routine
312	if no bins are given (bin_no, max_dist).
313
314	Returns
315	-------
316	bin_centers : (n), :class:`numpy.ndarray`
317	The bin centers.
318	gamma : (n) or (d, n), :class:`numpy.ndarray`
319	The estimated variogram values at bin centers.
320	Is stacked if multiple `directions` (d>1) are given.
321	counts : (n) or (d, n), :class:`numpy.ndarray`, optional
322	The number of point pairs found for each bin.
323	Is stacked if multiple `directions` (d>1) are given.
324	Only provided if `return_counts` is True.
325	normalizer : :any:`Normalizer`, optional
326	The fitted normalizer for the given data.
327	Only provided if `fit_normalizer` is True.
328
329	Notes
330	-----
331	Internally uses double precision and also returns doubles.
332
333	References
334	----------
335	.. [Webster2007] Webster, R. and Oliver, M. A.
336	"Geostatistics for environmental scientists.",
337	John Wiley & Sons. (2007)
338	"""
339	if bin_edges is not None:	1✔
340	bin_edges = np.atleast_1d(np.asarray(bin_edges, dtype=np.double))	1✔
341	bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0	1✔
342	# allow multiple fields at same positions (ndmin=2: first axis -> field ID)
343	# need to convert to ma.array, since list of ma.array is not recognised
344	field = np.ma.array(field, ndmin=2, dtype=np.double, copy=True)	1✔
345	masked = np.ma.is_masked(field) or np.any(mask)	1✔
346	# catch special case if everything is masked
347	if masked and np.all(mask):	1✔
348	bin_centers = np.empty(0) if bin_edges is None else bin_centers	1✔
349	estimates = np.zeros_like(bin_centers)	1✔
350	if return_counts:	1✔
351	return bin_centers, estimates, np.zeros_like(estimates, dtype=int)	1✔
352	return bin_centers, estimates	1✔
353	if not masked:	1✔
354	field = field.filled()	1✔
355	# check mesh shape
356	if mesh_type != "unstructured":	1✔
357	pos, __, dim = format_struct_pos_shape(	1✔
358	pos, field.shape, check_stacked_shape=True
359	)
360	pos = generate_grid(pos)	1✔
361	else:
362	pos, __, dim = format_unstruct_pos_shape(	1✔
363	pos, field.shape, check_stacked_shape=True
364	)
365	if latlon and dim != 2:	1✔
366	raise ValueError("Variogram: given field needs to be 2D for lat-lon.")	1✔
367	# prepare the field
368	pnt_cnt = len(pos[0])	1✔
369	field = field.reshape((-1, pnt_cnt))	1✔
370	# apply mask if wanted
371	if masked:	1✔
372	# if fields have different masks, take the minimal common mask
373	# given mask will be applied in addition
374	# selected region is the inverted masked (unmasked values)
375	if np.size(mask) > 1: # not only np.ma.nomask	1✔
376	select = np.invert(	1✔
377	np.logical_or(
378	np.reshape(mask, pnt_cnt), np.all(field.mask, axis=0)
379	)
380	)
381	else:
382	select = np.invert(np.all(field.mask, axis=0))	1✔
383	pos = pos[:, select]	1✔
384	field.fill_value = np.nan # use no-data val. for remaining masked vals	1✔
385	field = field[:, select].filled() # convert to ndarray	1✔
386	select = mask = None # free space	1✔
387	pnt_cnt = len(pos[0]) # pnt cnt reduced now	1✔
388	# set no_data values
389	if not np.isnan(no_data):	1✔
390	field[np.isclose(field, float(no_data))] = np.nan	1✔
391	# set directions
392	dir_no = 0	1✔
393	if direction is not None and dim > 1:	1✔
394	direction = np.atleast_2d(np.asarray(direction, dtype=np.double))	1✔
395	if len(direction.shape) > 2:	1✔
396	raise ValueError(f"Can't interpret directions: {direction}")	1✔
397	if direction.shape[1] != dim:	1✔
398	raise ValueError(f"Can't interpret directions: {direction}")	1✔
399	dir_no = direction.shape[0]	1✔
400	# convert given angles to direction vector
401	if angles is not None and direction is None and dim > 1:	1✔
402	direction = ang2dir(angles=angles, dtype=np.double, dim=dim)	1✔
403	dir_no = direction.shape[0]	1✔
404	# prepare directional variogram
405	if dir_no > 0:	1✔
406	if latlon:	1✔
407	raise ValueError("Directional variogram not allowed for lat-lon.")	1✔
408	norms = np.linalg.norm(direction, axis=1)	1✔
409	if np.any(np.isclose(norms, 0)):	1✔
410	raise ValueError(f"Zero length directions: {direction}")	1✔
411	# only unit-vectors for directions
412	direction = np.divide(direction, norms[:, np.newaxis])	1✔
413	# negative bandwidth to turn it off
414	bandwidth = float(bandwidth) if bandwidth is not None else -1.0	1✔
415	angles_tol = float(angles_tol)	1✔
416	# prepare sampled variogram
417	if sampling_size is not None and sampling_size < pnt_cnt:	1✔
418	sampled_idx = np.random.RandomState(sampling_seed).choice(	1✔
419	np.arange(pnt_cnt), sampling_size, replace=False
420	)
421	field = field[:, sampled_idx]	1✔
422	pos = pos[:, sampled_idx]	1✔
423	# create bins
424	if bin_edges is None:	1✔
425	bin_edges = standard_bins(	1✔
426	pos, dim, latlon, geo_scale=geo_scale, **std_bins
427	)
428	bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0	1✔
429	if latlon:	1✔
430	# internally we always use radians
431	bin_edges /= geo_scale	1✔
432	# normalize field
433	norm_field_out = remove_trend_norm_mean(	1✔
434	*(pos, field, mean, normalizer, trend),
435	check_shape=False,
436	stacked=True,
437	fit_normalizer=fit_normalizer,
438	)
439	field = norm_field_out[0] if fit_normalizer else norm_field_out	1✔
440	norm_out = (norm_field_out[1],) if fit_normalizer else ()	1✔
441	# select variogram estimator
442	cython_estimator = _set_estimator(estimator)	1✔
443	# run
444	if dir_no == 0:	1✔
445	# "h"aversine or "e"uclidean distance type
446	distance_type = "h" if latlon else "e"	1✔
447	estimates, counts = _unstructured(	1✔
448	field,
449	bin_edges,
450	pos,
451	estimator_type=cython_estimator,
452	distance_type=distance_type,
453	num_threads=config.NUM_THREADS,
454	)
455	else:
456	estimates, counts = _directional(	1✔
457	field,
458	bin_edges,
459	pos,
460	direction,
461	angles_tol,
462	bandwidth,
463	separate_dirs=_separate_dirs_test(direction, angles_tol),
464	estimator_type=cython_estimator,
465	num_threads=config.NUM_THREADS,
466	)
467	if dir_no == 1:	1✔
468	estimates, counts = estimates[0], counts[0]	1✔
469	est_out = (estimates, counts)	1✔
470	return (bin_centers,) + est_out[: 2 if return_counts else 1] + norm_out	1✔
471
472
473	def vario_estimate_axis(	1✔
474	field, direction="x", estimator="matheron", no_data=np.nan
475	):
476	r"""Estimates the variogram along array axis.
477
478	The indices of the given direction are used for the bins.
479	Uniform spacings along the given axis are assumed.
480
481	The algorithm calculates following equation:
482
483	.. math::
484	\gamma(r_k) = \frac{1}{2 N(r_k)} \sum_{i=1}^{N(r_k)} (z(\mathbf x_i) -
485	z(\mathbf x_i'))^2 \; ,
486
487	with :math:`r_k \leq \\| \mathbf x_i - \mathbf x_i' \\| < r_{k+1}`
488	being the bins.
489
490	Or if the estimator "cressie" was chosen:
491
492	.. math::
493	\gamma(r_k) = \frac{\frac{1}{2}\left(\frac{1}{N(r_k)}\sum_{i=1}^{N(r_k)}
494	\left\|z(\mathbf x_i) - z(\mathbf x_i')\right\|^{0.5}\right)^4}
495	{0.457 + 0.494 / N(r_k) + 0.045 / N^2(r_k)} \; ,
496
497	with :math:`r_k \leq \\| \mathbf x_i - \mathbf x_i' \\| < r_{k+1}`
498	being the bins.
499	The Cressie estimator is more robust to outliers [Webster2007]_.
500
501	Parameters
502	----------
503	field : :class:`numpy.ndarray` or :class:`numpy.ma.MaskedArray`
504	the spatially distributed data (can be masked)
505	direction : :class:`str` or :class:`int`
506	the axis over which the variogram will be estimated (x, y, z)
507	or (0, 1, 2, ...)
508	estimator : :class:`str`, optional
509	the estimator function, possible choices:
510
511	* "matheron": the standard method of moments of Matheron
512	* "cressie": an estimator more robust to outliers
513
514	Default: "matheron"
515
516	no_data : :class:`float`, optional
517	Value to identify missing data in the given field.
518	Default: `numpy.nan`
519
520	Returns
521	-------
522	:class:`numpy.ndarray`
523	the estimated variogram along the given direction.
524
525	Warnings
526	--------
527	It is assumed that the field is defined on an equidistant Cartesian grid.
528
529	Notes
530	-----
531	Internally uses double precision and also returns doubles.
532
533	References
534	----------
535	.. [Webster2007] Webster, R. and Oliver, M. A.
536	"Geostatistics for environmental scientists.",
537	John Wiley & Sons. (2007)
538	"""
539	missing_mask = (	1✔
540	np.isnan(field) if np.isnan(no_data) else np.isclose(field, no_data)
541	)
542	missing = np.any(missing_mask)	1✔
543	masked = np.ma.is_masked(field) or missing	1✔
544	if masked:	1✔
545	field = np.ma.array(field, ndmin=1, dtype=np.double)	1✔
546	if missing:	1✔
547	field.mask = np.logical_or(field.mask, missing_mask)	1✔
548	mask = np.ma.getmaskarray(field)	1✔
549	else:
550	field = np.atleast_1d(np.asarray(field, dtype=np.double))	1✔
551	missing_mask = None # free space	1✔
552
553	axis_to_swap = AXIS_DIR[direction] if direction in AXIS else int(direction)	1✔
554	# desired axis first, convert to 2D array afterwards
555	field = field.swapaxes(0, axis_to_swap)	1✔
556	field = field.reshape((field.shape[0], -1))	1✔
557	if masked:	1✔
558	mask = mask.swapaxes(0, axis_to_swap)	1✔
559	mask = mask.reshape((mask.shape[0], -1))	1✔
560
561	cython_estimator = _set_estimator(estimator)	1✔
562
563	if masked:	1✔
564	return _ma_structured(	1✔
565	field, mask, cython_estimator, num_threads=config.NUM_THREADS
566	)
567	return _structured(field, cython_estimator, num_threads=config.NUM_THREADS)	1✔
568
569
570	# for backward compatibility
571	vario_estimate_unstructured = vario_estimate	1✔
572	vario_estimate_structured = vario_estimate_axis	1✔

GeoStat-Framework / GSTools / 16463739022

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous