invert4geom.inversion

Module Contents

Functions

grav_column_der(grav_easting, grav_northing, ...)	Function to calculate the vertical derivate of the gravitational acceleration at
jacobian_annular(grav_easting, grav_northing, ...)	Function to calculate the Jacobian matrix using the annular cylinder
prism_properties(prisms_layer[, method])	extract prism properties from prism layer
jacobian_prism(prisms_properties, grav_easting, ...)	Function to calculate the Jacobian matrix with the vertical gravity derivative
jacobian(deriv_type, coordinates[, empty_jac, ...])	dispatcher for creating the jacobian matrix with 2 method options
solver(jac, residuals[, damping, solver_type])	Calculate shift to add to prism's for each iteration of the inversion. Finds
update_l2_norms(rmse, l2_norm)	update the l2 norm and delta l2 norm of the misfit
end_inversion(iteration_number, max_iterations, ...[, ...])	check if the inversion should be terminated
update_gravity_and_misfit(gravity_df, prisms_ds, ...)	calculate the forward gravity of the supplied prism layer, add the results to a
run_inversion(input_grav, input_grav_column, ...[, ...])	perform a geometric inversion, where the topography is updated to minimize the

grav_column_der(grav_easting, grav_northing, grav_upward, prism_easting, prism_northing, prism_top, prism_spacing, prism_density)

Function to calculate the vertical derivate of the gravitational acceleration at an observation point caused by a right, rectangular prism. Approximated with Hammer’s annulus approximation [1].

Parameters:

• grav_easting (NDArray) – coordinates of gravity observation points.

• grav_northing (NDArray) – coordinates of gravity observation points.

• grav_upward (NDArray) – coordinates of gravity observation points.

• prism_easting (NDArray) – coordinates of prism’s center in northing, easting, and upward directions, respectively

• prism_northing (NDArray) – coordinates of prism’s center in northing, easting, and upward directions, respectively

• prism_top (NDArray) – coordinates of prism’s center in northing, easting, and upward directions, respectively

• prism_spacing (float) – resolution of prism layer in meters

• prism_density (NDArray) – density of prisms, in kg/m^3

Returns:

array of vertical derivative of gravity at observation point for series of prisms

Return type:

NDArray

References

[1] (1,2)

Jack McCubbine. Airborne Gravity Across New Zealand - For An Improved Vertical Datum. PhD thesis, Victoria University of Wellington, 2016.

jacobian_annular(grav_easting, grav_northing, grav_upward, prism_easting, prism_northing, prism_top, prism_density, prism_spacing, jac)

Function to calculate the Jacobian matrix using the annular cylinder approximation. The resulting Jacobian is a matrix (numpy array) with a row per gravity observation and a column per prism. This approximates the prisms as an annulus [1], and calculates it’s vertical gravity derivative. Takes arrays from jacobian, feeds them into grav_column_der, and returns the jacobian.

Parameters:

• grav_easting (NDArray) – coordinates of gravity observation points

• grav_northing (NDArray) – coordinates of gravity observation points

• grav_upward (NDArray) – coordinates of gravity observation points

• prism_easting (NDArray) – coordinates of prism’s center in northing, easting, and upward directions, respectively

• prism_northing (NDArray) – coordinates of prism’s center in northing, easting, and upward directions, respectively

• prism_top (NDArray) – coordinates of prism’s center in northing, easting, and upward directions, respectively

• prism_density (NDArray) – density of prisms, in kg/m^3

• prism_spacing (float) – resolution of prism layer in meters

• jac (NDArray) – empty jacobian matrix with a row per gravity observation and a column per prism

Returns:

returns a jacobian matrix of shape (number of gravity points, number of prisms)

Return type:

NDArray

References

prism_properties(prisms_layer, method='itertools')

extract prism properties from prism layer

Parameters:

• prisms_layer (xr.Dataset) – harmonica prism layer

• method (str, optional) – choice of method to extract properties, by default “itertools”

Returns:

array of prism properties

Return type:

NDArray

jacobian_prism(prisms_properties, grav_easting, grav_northing, grav_upward, delta, jac)

Function to calculate the Jacobian matrix with the vertical gravity derivative as a numerical approximation with small prisms

Takes arrays from jacobian and calculates the jacobian.

Parameters:

• prisms_properties (NDArray) – array of prism properties of shape (number of prisms, 7) with the 7 entries for each prism being: west, east, south, north, bottom, top, density

• grav_easting (NDArray) – coordinates of gravity observation points.

• grav_northing (NDArray) – coordinates of gravity observation points.

• grav_upward (NDArray) – coordinates of gravity observation points.

• delta (float) – thickness in meters of small prisms used to calculate vertical derivative

• jac (NDArray) – empty jacobian matrix with a row per gravity observation and a column per prism

Returns:

returns a NDArray of shape (number of gravity points, number of prisms)

Return type:

NDArray

jacobian(deriv_type, coordinates, empty_jac=None, prisms_layer=None, prism_spacing=None, prism_size=None, prisms_properties_method='itertools')

dispatcher for creating the jacobian matrix with 2 method options

Parameters:

• deriv_type (str) – choose between “annulus” and “prisms” methods of calculating the vertical derivative of gravity of a prism

• coordinates (pd.DataFrame) – coordinate dataframe of gravity observation points with columns “easting”, “northing”, “upward”

• empty_jac (NDArray, optional) – optionally provide an empty jacobian matrix of shape (number of gravity observations x number of prisms), by default None

• prisms_layer (xr.Dataset, optional) – harmonica prism layer, by default None

• prism_spacing (float, optional) – resolution of prism layer, by default None

• prism_size (float, optional) – height of prisms for small prism vertical derivative method, by default None

• prisms_properties_method (str, optional) – method for extracting prism properties, by default “itertools”

Returns:

a filled out jacobian matrix

Return type:

NDArray

solver(jac, residuals, damping=None, solver_type='scipy least squares')

Calculate shift to add to prism’s for each iteration of the inversion. Finds the least-squares solution to the Jacobian and the gravity residual

Parameters:

• jac (NDArray) – input jacobian matrix with a row per gravity observation, and a column per prisms.

• residuals (NDArray) – array of gravity residuals

• damping (float | None, optional) – positive damping (Tikhonov 0th order) regularization

• solver_type – ‘verde least squares’, ‘scipy least squares’, ‘scipy conjugate’, ‘numpy least squares’, ‘steepest descent’, ‘gauss newton’, } optional choose which solving method to use, by default “scipy least squares”

Returns:

array of correction values to apply to each prism.

Return type:

NDArray

update_l2_norms(rmse, l2_norm)

update the l2 norm and delta l2 norm of the misfit

Parameters:

• rmse (float) – root mean square error of the residual gravity misfit

• l2_norm (float) – l2 norm of the residual gravity misfit

Returns:

updated l2 norm and delta l2 norm

Return type:

tuple[float, float]

end_inversion(iteration_number, max_iterations, l2_norm, starting_l2_norm, l2_norm_tolerance, delta_l2_norm, previous_delta_l2_norm, delta_l2_norm_tolerance, perc_increase_limit=0.2)

check if the inversion should be terminated

Parameters:

• iteration_number (int) – the iteration number, starting at 1 not 0

• max_iterations (int) – the maximum allowed iterations, inclusive and starting at 1

• l2_norm (float) – the current iteration’s l2 norm

• starting_l2_norm (float) – the l2 norm of iteration 1

• l2_norm_tolerance (float) – the l2 norm value to end the inversion at

• delta_l2_norm (float) – the current iteration’s delta l2 norm

• previous_delta_l2_norm (float) – the delta l2 norm of the previous iteration

• delta_l2_norm_tolerance (float) – the delta l2 norm value to end the inversion at

• perc_increase_limit (float, optional) – the set tolerance for decimal percentage increase relative to the starting l2 norm, by default 0.20

Returns:

first term is a boolean of whether or not to end the inversion, second term is a list of termination reasons.

Return type:

tuple[bool, list[str]]

update_gravity_and_misfit(gravity_df, prisms_ds, input_grav_column, iteration_number)

calculate the forward gravity of the supplied prism layer, add the results to a new dataframe column, and update the residual misfit. The supplied gravity dataframe needs a ‘reg’ column, which describes the regional component and can be 0.

Parameters:

• gravity_df (pd.DataFrame) – gravity dataframe with gravity observation coordinate columns (‘easting’, ‘northing’, ‘upwards’), a gravity data column, set by input_grav_column, and a regional gravity column (‘reg’).

• prisms_ds (xr.Dataset) – harmonica prism layer

• input_grav_column (str) – name of gravity data column

• iteration_number (int) – iteration number to use in updated column names

Returns:

a gravity dataframe with 2 new columns, one for the iterations forward gravity and one for the iterations residual misfit.

Return type:

pd.DataFrame

run_inversion(input_grav, input_grav_column, prism_layer, density_contrast, zref, max_iterations, l2_norm_tolerance=0.2, delta_l2_norm_tolerance=1.001, perc_increase_limit=0.1, deriv_type='annulus', jacobian_prism_size=1, solver_type='scipy least squares', solver_damping=None, upper_confining_layer=None, lower_confining_layer=None, weights_after_solving=False, inversion_region=None, plot_convergence=False, plot_dynamic_convergence=False)

perform a geometric inversion, where the topography is updated to minimize the residual misfit between the forward gravity of the layer, and the observed gravity. To aid in regularizing an ill-posed problem choose any of the following options: * add damping to the solver, with solver_damping * weight the surface correction values with a weighting grid with weights_after_solving and the weights variable of the prisms dataset * bound the topography of the layer, with upper_confining_layer and lower_confining_layer

Parameters:

• input_grav (pd.DataFrame) – dataframe with gravity data and coordinates, must have columns “res” and “reg” for residual and regional gravity, and coordinate columns “easting”, “northing”, and “upward”.

• input_grav_column (str) – column name containing the gravity data before regional separation

• prism_layer (xr.Dataset) – starting prism layer

• density_contrast (float) – density contrast of the prisms layer, should be same value used to create the starting model

• zref (float) – reference height of the prisms layer, should be same value used to create the starting model

• max_iterations (int) – the maximum allowed iterations, inclusive and starting at 1

• l2_norm_tolerance (float, optional) – _the l2 norm value to end the inversion at, by default 0.2

• delta_l2_norm_tolerance (float, optional) – the delta l2 norm value to end the inversion at, by default 1.001

• perc_increase_limit (float, optional) – the set tolerance for decimal percentage increase relative to the starting l2 norm, by default 0.10

• deriv_type (str, optional) – either “annulus” or “prism” to determine method of calculating the vertical derivate of gravity of a prism, by default “annulus”

• jacobian_prism_size (float, optional) – height of prisms in meters for vertical derivative, by default 1

• solver_type (str, optional) – solver type to use, by default “scipy least squares”

• solver_damping (float | None, optional) – damping parameter for regularization of the solver, by default None

• upper_confining_layer (xr.DataArray | None, optional) – topographic layer to use as upper limit for inverted topography, by default None

• lower_confining_layer (xr.DataArray | None, optional) – topographic layer to use as lower limit for inverted topography, by default None

• weights_after_solving (bool, optional) – use “weights” variable of prisms dataset to scale surface corrections grid, by default False, by default False

• inversion_region (tuple[float, float, float, float]) – inside region to calculated residual RMSE within, in the form (min_easting, max_easting, min_northing, max_northing)

• plot_convergence (bool, optional) – plot the misfit convergence, by default False

• plot_dynamic_convergence (bool, optional) – plot the misfit convergence dynamically, by default False

Returns:

prisms_df: pd.DataFrame, prism properties for each iteration, gravity: pd.DataFrame, gravity anomalies for each iteration, params: dict, Properties of the inversion such as kwarg values, elapsed_time: float, time in seconds for the inversion to run

Return type:

tuple[pd.DataFrame, pd.DataFrame, dict[str, Any], float]