localization

Localization tools, including distance and tapering comps.

A good introduction to localization: Sakov (2011), Computational Geosciences: 'Relation between two common localisation methods for the EnKF'.

dist2coeff(dists, radius, tag=None)

Compute tapering coefficients corresponding to a distances.

NB: The radius is internally adjusted such that, independently of 'tag', coeff==np.exp(-0.5) when distance==radius.

This is largely based on Sakov's enkf-matlab code. Two bugs have here been fixed: - The constants were slightly wrong, as noted in comments below. - It forgot to take sqrt() of coeffs when applying them through 'local analysis'.

inds_and_coeffs(dists, radius, cutoff=0.001, tag=None)

Compute indices and coefficients of localization.

• inds : the indices of pts that are "close to" centre.
• coeffs : the corresponding tapering coefficients.

nd_Id_localization(shape, batch_shape=None, obs_inds=None, periodic=True)

Localize Id (direct) point obs of an N-D, homogeneous, rectangular domain.

pairwise_distances(A, B=None, domain=None)

Euclidian distance (not squared) between pts. in A and B.

Parameters:

Name	Type	Description	Default
A	array of shape `(nPoints, nDims)`.	A collection of points.	required
B		Same as A, but nPoints can differ.	None
domain	tuple	Assume the domain is a periodic hyper-rectangle whose edges along dimension i span from 0 to domain[i]. NB: Behaviour not defined if any(A.max(0) > domain), and likewise for B.	None

Returns:

Type	Description
Array of of shape `(nPointsA, nPointsB)`.

Examples:

>>> A = [[0, 0], [0, 1], [1, 0], [1, 1]]
>>> with np.printoptions(precision=2):
...     print(pairwise_distances(A))
[[0.   1.   1.   1.41]
 [1.   0.   1.41 1.  ]
 [1.   1.41 0.   1.  ]
 [1.41 1.   1.   0.  ]]

The function matches pdist(..., metric='euclidean'), but is faster:

>>> from scipy.spatial.distance import pdist, squareform
>>> (pairwise_distances(A) == squareform(pdist(A))).all()
True

As opposed to pdist, it also allows comparing A to a different set of points, B, without the augmentation/block tricks needed for pdist.

>>> A = np.arange(4)[:, None]
>>> pairwise_distances(A, [[2]]).T
array([[2., 1., 0., 1.]])

Illustration of periodicity:

>>> pairwise_distances(A, domain=(4, ))
array([[0., 1., 2., 1.],
       [1., 0., 1., 2.],
       [2., 1., 0., 1.],
       [1., 2., 1., 0.]])

NB: If an input array is 1-dim, it is seen as a single point.

>>> pairwise_distances(np.arange(4))
array([[0.]])

rectangular_partitioning(shape, steps, do_ind=True)

N-D rectangular batch generation.

Parameters:

Name	Type	Description	Default
shape	len(grid[dim]) for dim in range(ndim)		required
steps	step_len[dim] for dim in range(ndim)		required

Returns:

Type	Description
A list of batches,
where each element (batch) is a list of indices.

Example

shape = [4, 13] ... batches = rectangular_partitioning(shape, [2, 4], do_ind=False) ... nB = len(batches) ... values = np.random.choice(np.arange(nB), nB, 0) ... Z = np.zeros(shape) ... for ib, b in enumerate(batches): ... Z[tuple(b)] = values[ib] ... plt.imshow(Z) # doctest: +SKIP