geoist.gridder package

Submodules

geoist.gridder.genpnt module

Name : genpnt.py Created on : 2018/11/24 08:57 Author : Steve Chen <chenshi@cea-igp.ac.cn> Affiliation : Institute of Geophysics, CEA. Version : 0.1.0 Copyright : Copyright (C) 2018-2020 GEOIST Development Team. All Rights Reserved. License : Distributed under the MIT License. See LICENSE.txt for more info. Github : https://igp-gravity.github.io/ Description : Generate points on a map as regular grids or points scatters.

geoist.gridder.genpnt.circular_scatter(area, n, z=None, random=False, seed=None)

Generate a set of n points positioned in a circular array.

The diameter of the circle is equal to the smallest dimension of the area

Parameters:

  • area : list = [x1, x2, y1, y2]
    Area inside of which the points are contained
  • n : int
    Number of points
  • z : float or 1d-array
    Optional. z coordinate of the points. If given, will return an array with the value z.
  • random : True or False
    If True, positions of the points on the circle will be chosen at random
  • seed : None or int
    Seed used to generate the pseudo-random numbers if random==True. If None, will use a different seed every time. Use the same seed to generate the same random sequence.

Returns:

  • [x, y]
    Numpy arrays with the x and y coordinates of the points
  • [x, y, z]
    If z given. Arrays with the x, y, and z coordinates of the points
geoist.gridder.genpnt.regular(area, shape, z=None)

Create a regular grid.

The x directions is North-South and y East-West. Imagine the grid as a matrix with x varying in the lines and y in columns.

Returned arrays will be flattened to 1D with numpy.ravel.

Warning

As of version 0.4, the shape argument was corrected to be shape = (nx, ny) instead of shape = (ny, nx).

Parameters:

  • area
    (x1, x2, y1, y2): Borders of the grid
  • shape
    Shape of the regular grid, ie (nx, ny).
  • z
    Optional. z coordinate of the grid points. If given, will return an array with the value z.

Returns:

  • [x, y]
    Numpy arrays with the x and y coordinates of the grid points
  • [x, y, z]
    If z given. Numpy arrays with the x, y, and z coordinates of the grid points

Examples:

>>> x, y = regular((0, 10, 0, 5), (5, 3))
>>> print(x)
[  0.    0.    0.    2.5   2.5   2.5   5.    5.    5.    7.5   7.5   7.5 10.   10.   10. ]
>>> print(x.reshape((5, 3)))
[[  0.    0.    0. ]
 [  2.5   2.5   2.5]
 [  5.    5.    5. ]
 [  7.5   7.5   7.5]
 [ 10.   10.   10. ]]
>>> print(y.reshape((5, 3)))
[[ 0.   2.5  5. ]
 [ 0.   2.5  5. ]
 [ 0.   2.5  5. ]
 [ 0.   2.5  5. ]
 [ 0.   2.5  5. ]]
>>> x, y = regular((0, 0, 0, 5), (1, 3))
>>> print(x.reshape((1, 3)))
[[ 0.  0.  0.]]
>>> print(y.reshape((1, 3)))
[[ 0.   2.5  5. ]]
>>> x, y, z = regular((0, 10, 0, 5), (5, 3), z=-10)
>>> print(z.reshape((5, 3)))
[[-10. -10. -10.]
 [-10. -10. -10.]
 [-10. -10. -10.]
 [-10. -10. -10.]
 [-10. -10. -10.]]
geoist.gridder.genpnt.scatter(area, n, z=None, seed=None)

Create an irregular grid with a random scattering of points.

Parameters:

  • area
    (x1, x2, y1, y2): Borders of the grid
  • n
    Number of points
  • z
    Optional. z coordinate of the points. If given, will return an array with the value z.
  • seed : None or int
    Seed used to generate the pseudo-random numbers. If None, will use a different seed every time. Use the same seed to generate the same random points.

Returns:

  • [x, y]
    Numpy arrays with the x and y coordinates of the points
  • [x, y, z]
    If z given. Arrays with the x, y, and z coordinates of the points

Examples:

>>> # Passing in a seed value will ensure that scatter will return the same
>>> # values given the same input. Use seed=None if you don't want this.
>>> x, y = scatter((0, 10, 0, 2), 4, seed=0)
>>> # Small function to print the array values with 4 decimal places
>>> pprint = lambda arr: print(', '.join('{:.4f}'.format(i) for i in arr))
>>> pprint(x)
5.4881, 7.1519, 6.0276, 5.4488
>>> pprint(y)
0.8473, 1.2918, 0.8752, 1.7835

>>> # scatter can take the z argument as well
>>> x2, y2, z2 = scatter((-10, 1, 1245, 3556), 6, z=-150, seed=2)
>>> pprint(x2)
-5.2041, -9.7148, -3.9537, -5.2115, -5.3760, -6.3663
>>> pprint(y2)
1717.9430, 2676.1352, 1937.5020, 1861.6378, 2680.4403, 2467.8474
>>> pprint(z2)
-150.0000, -150.0000, -150.0000, -150.0000, -150.0000, -150.0000

geoist.gridder.interpolation module

Name : interpolation.py Created on : 2018/11/24 08:57 Author : Steve Chen <chenshi@cea-igp.ac.cn> Affiliation : Institute of Geophysics, CEA. Version : 0.1.0 Copyright : Copyright (C) 2018-2020 GEOIST Development Team. All Rights Reserved. License : Distributed under the MIT License. See LICENSE.txt for more info. Github : https://igp-gravity.github.io/ Description : 2D interpolation, griding, and profile extraction.

geoist.gridder.interpolation.fill_nans(x, y, v, xp, yp, vp)

” Fill in the NaNs or masked values on interpolated points using nearest neighbors.

Warning

Operation is performed in place. Replaces the NaN or masked values of the original array!

Parameters:

  • x, y : 1D arrays
    Arrays with the x and y coordinates of the original data points (not interpolated).
  • v : 1D array
    Array with the scalar value assigned to the data points (not interpolated).
  • xp, yp : 1D arrays
    Points where the data values were interpolated.
  • vp : 1D array
    Interpolated data values (the one that has NaNs or masked values to replace).
geoist.gridder.interpolation.interp(x, y, v, shape, area=None, algorithm='cubic', extrapolate=False)

Interpolate spacial data onto a regular grid.

Utility function that generates a regular grid with regular() and calls interp_at() on the generated points.

Parameters:

  • x, y : 1D arrays
    Arrays with the x and y coordinates of the data points.
  • v : 1D array
    Array with the scalar value assigned to the data points.
  • shape : tuple = (nx, ny)
    Shape of the interpolated regular grid, ie (nx, ny).
  • area : tuple = (x1, x2, y1, y2)
    The are where the data will be interpolated. If None, then will get the area from x and y.
  • algorithm : string
    Interpolation algorithm. Either 'cubic', 'nearest', 'linear' (see scipy.interpolate.griddata).
  • extrapolate : True or False
    If True, will extrapolate values outside of the convex hull of the data points.

Returns:

  • [x, y, v]
    Three 1D arrays with the interpolated x, y, and v
geoist.gridder.interpolation.interp_at(x, y, v, xp, yp, algorithm='cubic', extrapolate=False)

Interpolate spacial data onto specified points.

Wraps scipy.interpolate.griddata.

Parameters:

  • x, y : 1D arrays
    Arrays with the x and y coordinates of the data points.
  • v : 1D array
    Array with the scalar value assigned to the data points.
  • xp, yp : 1D arrays
    Points where the data values will be interpolated
  • algorithm : string
    Interpolation algorithm. Either 'cubic', 'nearest', 'linear' (see scipy.interpolate.griddata)
  • extrapolate : True or False
    If True, will extrapolate values outside of the convex hull of the data points.

Returns:

  • v : 1D array
    1D array with the interpolated v values.
geoist.gridder.interpolation.profile(x, y, v, point1, point2, size, algorithm='cubic')

Extract a profile between 2 points from spacial data.

Uses interpolation to calculate the data values at the profile points.

Parameters:

  • x, y : 1D arrays
    Arrays with the x and y coordinates of the data points.
  • v : 1D array
    Array with the scalar value assigned to the data points.
  • point1, point2 : lists = [x, y]
    Lists the x, y coordinates of the 2 points between which the profile will be extracted.
  • size : int
    Number of points along the profile.
  • algorithm : string
    Interpolation algorithm. Either 'cubic', 'nearest', 'linear' (see scipy.interpolate.griddata).

Returns:

  • [xp, yp, distances, vp] : 1d arrays
    xp and yp are the x, y coordinates of the points along the profile. distances are the distances of the profile points from point1. vp are the data points along the profile.

geoist.gridder.padding module

Name : padding.py Created on : 2018/11/24 08:57 Author : Steve Chen <chenshi@cea-igp.ac.cn> Affiliation : Institute of Geophysics, CEA. Version : 0.1.0 Copyright : Copyright (C) 2018-2020 GEOIST Development Team. All Rights Reserved. License : Distributed under the MIT License. See LICENSE.txt for more info. Github : https://igp-gravity.github.io/ Description : Apply padding to data grids using different algorithms for the filling.

geoist.gridder.padding.pad_array(a, npd=None, padtype='OddReflectionTaper')

Return a padded array of arbitrary dimension.

The function takes an array of arbitrary dimension and pads it either to the dimensions given by the tuple npd, or to the next power of 2 if npd is not given.

An odd reflection with a cosine taper (padtype='OddReflectionTaper') is the preferred method of padding for Fourier Transform operations. The odd reflection optimally preserves the frequency content while adding minimal sharp inflections.

Note

Requires gridded data of the same dimension as the desired output (i.e. no flattened arrays; use reshape).

Note

This function returns a deep copy of the original array.

Parameters:

  • a : array
    Array (N-D) to be padded
  • npd : tuple (optional)
    Desired shape of new padded array. If not provided, the nearest power of 2 will be used.
  • padtype : string (optional)
    What method will be used to pad the new values. Can be lower or upper case. Options:
    • oddreflectiontaper: Generates odd reflection then tapers to the mean using a cosine function (Default)
    • oddreflection: Pads with the odd reflection, with no taper
    • reflection: Pads with simple reflection
    • lintaper: Linearly tapers to the mean
    • value: Numeric value (e.g., '10.4'). Input a float or integer directly.
    • edge: Uses the edge value as a constant pad
    • mean: Uses the mean of the vector along each axis

Returns:

  • ap : numpy array
    Padded array. The array core is a deep copy of the original array
  • nps : list
    List of tuples containing the number of elements padded onto each dimension.

Examples:

>>> import numpy as np
>>> z = np.array([3, 4, 4, 5, 6])
>>> zpad, nps = pad_array(z)
>>> print(zpad)
[ 4.4  3.2  3.   4.   4.   5.   6.   4.4]
>>> print(nps)
[(2, 1)]

>>> shape = (5, 6)
>>> z = np.ones(shape)
>>> zpad, nps = pad_array(z, padtype='5')
>>> zpad
array([[ 5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.],
       [ 5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.],
       [ 5.,  1.,  1.,  1.,  1.,  1.,  1.,  5.],
       [ 5.,  1.,  1.,  1.,  1.,  1.,  1.,  5.],
       [ 5.,  1.,  1.,  1.,  1.,  1.,  1.,  5.],
       [ 5.,  1.,  1.,  1.,  1.,  1.,  1.,  5.],
       [ 5.,  1.,  1.,  1.,  1.,  1.,  1.,  5.],
       [ 5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.]])
>>> print(nps)
[(2, 1), (1, 1)]
geoist.gridder.padding.pad_coords(xy, shape, nps)

Apply padding to coordinate vectors.

Designed to be used in concert with pad_array(), this function takes a list of coordinate vectors and pads them using the same number of elements as the padding of the data array.

Note

This function returns a list of arrays in the same format as, for example, regular(). It is a list of flattened np.meshgrid for each vector in the same order as was input through argument xy.

Parameters:

  • xy : list
    List of arrays of coordinates
  • shape : tuple
    Size of original array
  • nps : list
    List of tuples containing the number of elements padded onto each dimension (use the output from pad_array()).

Returns:

  • coordspad : list
    List of padded coordinate arrays

Examples:

>>> import numpy as np
>>> from fatiando.gridder import regular
>>> shape = (5, 6)
>>> x, y, z = regular((-10, 10, -20, 0), shape, z=-25)
>>> gz = np.zeros(shape)
>>> gzpad, nps = pad_array(gz)
>>> print(x.reshape(shape)[:, 0])
[-10.  -5.   0.   5.  10.]
>>> print(y.reshape(shape)[0, :])
[-20. -16. -12.  -8.  -4.   0.]
>>> xy = [x, y]
>>> N = pad_coords(xy, shape, nps)
>>> print(N[0].reshape(gzpad.shape)[:, 0])
[-20. -15. -10.  -5.   0.   5.  10.  15.]
>>> print(N[1].reshape(gzpad.shape)[0, :])
[-24. -20. -16. -12.  -8.  -4.   0.   4.]
geoist.gridder.padding.unpad_array(a, nps)

Remove padding from an array.

This function takes a padded array and removes the padding from both sides. Designed to use the output of pad_array().

Note

Unlike pad_array(), this function returns a slice of the input array. Therefore, any changes to the padded array will be reflected in the unpadded array!

Parameters:

  • a : array
    Array to be un-padded. Can be of arbitrary dimension.
  • nps : list
    List of tuples giving the min and max indices for the cutoff. Use the value returned by pad_array().

Returns:

  • b : array
    Array of same dimension as a, with padding removed

Examples:

>>> import numpy as np
>>> z = np.array([3, 4, 4, 5, 6])
>>> zpad, nps = pad_array(z)
>>> print(zpad)
[ 4.4  3.2  3.   4.   4.   5.   6.   4.4]
>>> zunpad = unpad_array(zpad, nps)
>>> print(zunpad)
[ 3.  4.  4.  5.  6.]

geoist.gridder.slicing module

Name : slicing.py Created on : 2018/11/24 08:57 Author : Steve Chen <chenshi@cea-igp.ac.cn> Affiliation : Institute of Geophysics, CEA. Version : 0.1.0 Copyright : Copyright (C) 2018-2020 GEOIST Development Team. All Rights Reserved. License : Distributed under the MIT License. See LICENSE.txt for more info. Github : https://igp-gravity.github.io/ Description : Functions for segmenting spacial data (windowing, cutting, etc).

geoist.gridder.slicing.cut(x, y, scalars, area)

Return a subsection of a grid.

The returned subsection is not a copy! In technical terms, returns a slice of the numpy arrays. So changes made to the subsection reflect on the original grid. Use numpy.copy to make copies of the subsections and avoid this.

Parameters:

  • x, y
    Arrays with the x and y coordinates of the data points.
  • scalars
    List of arrays with the scalar values assigned to the grid points.
  • area
    (x1, x2, y1, y2): Borders of the subsection

Returns:

  • [subx, suby, subscalars]
    Arrays with x and y coordinates and scalar values of the subsection.

Examples:

>>> import numpy as np
>>> x = np.array([1, 2, 3, 4, 5, 6])
>>> y = np.array([10, 11, 12, 13, 14, 15])
>>> data = np.array([42, 65, 92, 24, 135, 456])
>>> area = [2.5, 5.5, 12, 15]
>>> xs, ys, [datas] = cut(x, y, [data], area)
>>> print(xs)
[3 4 5]
>>> print(ys)
[12 13 14]
>>> print(datas)
[ 92  24 135]
>>> # This also works for 2D-arrays
>>> x = np.array([[1, 1, 1],
...               [2, 2, 2],
...               [3, 3, 3]])
>>> y = np.array([[5, 7, 9],
...               [5, 7, 9],
...               [5, 7, 9]])
>>> data = np.array([[12, 84, 53],
...                  [43, 79, 29],
...                  [45, 27, 10]])
>>> area = [0.5, 2.5, 7, 9]
>>> xs, ys, [datas] = cut(x, y, [data], area)
>>> print(xs)
[1 1 2 2]
>>> print(ys)
[7 9 7 9]
>>> print(datas)
[84 53 79 29]
geoist.gridder.slicing.inside(x, y, area)

Tell which indices of an array fall inside an area.

Parameters:

  • x, y : ndarrays
    The x and y coordinate vectors.
  • area : list = [xmin, xmax, ymin, ymax]
    x and y limits of the area.

Returns:

  • is_inside : ndarray of booleans
    An array of booleans. Will be True if the respective coordinates fall inside the area, False otherwise.

Examples:

>>> import numpy as np
>>> x = np.array([1, 2, 3, 4, 5, 6])
>>> y = np.array([10, 11, 12, 13, 14, 15])
>>> area = [2.5, 5.5, 12, 15]
>>> is_inside = inside(x, y, area)
>>> print(is_inside)
[False False  True  True  True False]
>>> # This also works for 2D-arrays
>>> x = np.array([[1, 1, 1],
...               [2, 2, 2],
...               [3, 3, 3]])
>>> y = np.array([[5, 7, 9],
...               [5, 7, 9],
...               [5, 7, 9]])
>>> area = [0.5, 2.5, 7, 9]
>>> is_inside = inside(x, y, area)
>>> print(is_inside)
[[False  True  True]
 [False  True  True]
 [False False False]]

Module contents

Create and operate on data grids, scatters, and profiles.