invert4geom.inversion#
Functions#
|
Function to calculate the vertical derivate of the gravitational acceleration at |
|
Function to calculate the Jacobian matrix using the annular cylinder |
|
extract prism properties from prism layer |
|
Function to calculate the Jacobian matrix with the vertical gravity derivative |
|
dispatcher for creating the jacobian matrix with 2 method options |
|
Calculate shift to add to prism's for each iteration of the inversion. Finds |
|
update the l2 norm and delta l2 norm of the misfit. |
|
check if the inversion should be terminated |
|
calculate the forward gravity of the supplied prism layer, add the results to a |
|
perform a geometric inversion, where the topography is updated to minimize the |
|
This function runs the full inversion workflow. Depending on the input parameters, |
Module Contents#
- grav_column_der(grav_easting, grav_northing, grav_upward, prism_easting, prism_northing, prism_top, prism_spacing, prism_density)[source]#
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 (numpy.ndarray) – coordinates of gravity observation points.
grav_northing (numpy.ndarray) – coordinates of gravity observation points.
grav_upward (numpy.ndarray) – coordinates of gravity observation points.
prism_easting (numpy.ndarray) – coordinates of prism’s center in northing, easting, and upward directions, respectively
prism_northing (numpy.ndarray) – coordinates of prism’s center in northing, easting, and upward directions, respectively
prism_top (numpy.ndarray) – coordinates of prism’s center in northing, easting, and upward directions, respectively
prism_spacing (float) – resolution of prism layer in meters
prism_density (numpy.ndarray) – density of prisms, in kg/m^3
- Returns:
array of vertical derivative of gravity at observation point for series of prisms
- Return type:
References
- jacobian_annular(grav_easting, grav_northing, grav_upward, prism_easting, prism_northing, prism_top, prism_density, prism_spacing, jac)[source]#
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 (numpy.ndarray) – coordinates of gravity observation points
grav_northing (numpy.ndarray) – coordinates of gravity observation points
grav_upward (numpy.ndarray) – coordinates of gravity observation points
prism_easting (numpy.ndarray) – coordinates of prism’s center in northing, easting, and upward directions, respectively
prism_northing (numpy.ndarray) – coordinates of prism’s center in northing, easting, and upward directions, respectively
prism_top (numpy.ndarray) – coordinates of prism’s center in northing, easting, and upward directions, respectively
prism_density (numpy.ndarray) – density of prisms, in kg/m^3
prism_spacing (float) – resolution of prism layer in meters
jac (numpy.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:
References
- _prism_properties(prisms_layer, method='itertools')[source]#
extract prism properties from prism layer
- Parameters:
prisms_layer (xarray.Dataset) – harmonica prism layer
method (str, optional) – choice of method to extract properties, by default “itertools”
- Returns:
array of prism properties
- Return type:
- jacobian_prism(prisms_properties, grav_easting, grav_northing, grav_upward, delta, jac)[source]#
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 (numpy.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 (numpy.ndarray) – coordinates of gravity observation points.
grav_northing (numpy.ndarray) – coordinates of gravity observation points.
grav_upward (numpy.ndarray) – coordinates of gravity observation points.
delta (float) – thickness in meters of small prisms used to calculate vertical derivative
jac (numpy.ndarray) – empty jacobian matrix with a row per gravity observation and a column per prism
- Returns:
returns a numpy.ndarray of shape (number of gravity points, number of prisms)
- Return type:
- jacobian(deriv_type, coordinates, empty_jac=None, prisms_layer=None, prism_spacing=None, prism_size=None, prisms_properties_method='itertools')[source]#
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 (pandas.DataFrame) – coordinate dataframe of gravity observation points with columns “easting”, “northing”, “upward”
empty_jac (numpy.ndarray, optional) – optionally provide an empty jacobian matrix of shape (number of gravity observations x number of prisms), by default None
prisms_layer (xarray.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:
- solver(jac, residuals, damping=None, solver_type='scipy least squares')[source]#
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 (numpy.ndarray) – input jacobian matrix with a row per gravity observation, and a column per prisms.
residuals (numpy.ndarray) – array of gravity residuals
damping (float | None, optional) – positive damping (Tikhonov 0th order) regularization
solver_type (str, optional) – choose which solving method to use, by default “scipy least squares”
- Returns:
array of correction values to apply to each prism.
- Return type:
- update_l2_norms(current_rmse, last_l2_norm)[source]#
update the l2 norm and delta l2 norm of the misfit.
We want the misfit (L2-norm) to be steadily decreasing with each iteration. If it increases, something is wrong, stop inversion. If it doesn’t decrease enough, inversion has finished and can be stopped. Delta L2 norm starts at +inf, and should decreases with each iteration. If it gets close to 1, the iterations aren’t making progress and can be stopped. A value of 1.001 means the L2 norm has only decrease by 0.1% between iterations and RMSE has only decreased by 0.05%.
- end_inversion(iteration_number, max_iterations, l2_norms, l2_norm_tolerance, delta_l2_norm, previous_delta_l2_norm, delta_l2_norm_tolerance, perc_increase_limit)[source]#
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_norms (float) – a list of each iteration’s l2 norm
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) – the set tolerance for decimal percentage increase relative to the starting l2 norm
- Returns:
end (bool) – whether or not to end the inversion
termination_reason (list[str]) – a list of termination reasons
- Return type:
- update_gravity_and_misfit(gravity_df, prisms_ds, iteration_number)[source]#
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 (pandas.DataFrame) – gravity dataframe with gravity observation coordinate columns (‘easting’, ‘northing’), a gravity data column ‘gravity_anomaly’, and a regional gravity column (‘reg’).
prisms_ds (xarray.Dataset) – harmonica prism layer
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:
- run_inversion(grav_df, prism_layer, max_iterations, l2_norm_tolerance=0.2, delta_l2_norm_tolerance=1.001, perc_increase_limit=0.2, deriv_type='annulus', jacobian_prism_size=1, solver_type='scipy least squares', solver_damping=None, upper_confining_layer=None, lower_confining_layer=None, apply_weighting_grid=False, weighting_grid=None, plot_convergence=False, plot_dynamic_convergence=False, results_fname=None, progressbar=True)[source]#
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 apply_weighting_grid and the weighting_grid argument * bound the topography of the layer, with upper_confining_layer and lower_confining_layer
- Parameters:
grav_df (pandas.DataFrame) – dataframe with gravity data and coordinates, must have columns “gravity_anomaly”, “misfit”, “res” and “reg”, and coordinate columns “easting” and “northing”.
prism_layer (xarray.Dataset) – starting prism layer
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 (xarray.DataArray | None, optional) – topographic layer to use as upper limit for inverted topography, by default None
lower_confining_layer (xarray.DataArray | None, optional) – topographic layer to use as lower limit for inverted topography, by default None
apply_weighting_grid (bool, optional) – use “weighting_grid” to scale surface corrections grid, by default False, by default False
plot_convergence (bool, optional) – plot the misfit convergence, by default False
plot_dynamic_convergence (bool, optional) – plot the misfit convergence dynamically, by default False
results_fname (str, optional) – filename to save results to, by default None
progressbar (bool, optional) – set to False to turn off the inversion progressbar
weighting_grid (xarray.DataArray | None)
- Returns:
prisms_df (pandas.DataFrame) – prism properties for each iteration, including columns starting_topo for the unchanged starting topography model and topo for the updated topography model
gravity (pandas.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[pandas.DataFrame, pandas.DataFrame, dict[str, Any], float]
- run_inversion_workflow(grav_df, create_starting_topography=False, create_starting_prisms=False, calculate_starting_gravity=False, calculate_regional_misfit=False, run_damping_cv=False, run_zref_or_density_cv=False, run_kfolds_zref_or_density_cv=False, plot_cv=False, fname=None, **kwargs)[source]#
This function runs the full inversion workflow. Depending on the input parameters, it will: 1) create a starting topography model 2) create a starting prism model 3) calculate the starting gravity of the prism model 4) calculate the gravity misfit 5) calculate the regional and residual components of the misfit 6) run the inversion you can choose to run cross validations for damping, density, and zref
- Parameters:
grav_df (pandas.DataFrame) – gravity dataframe with gravity data, must have coordinate columns ‘easting’, and ‘northing’. It must also have a gravity data column ‘gravity_anomaly’. Optionally should have columns ‘starting_gravity’, ‘misfit’, ‘reg’, ‘res’.
create_starting_topography (bool, optional) – Choose whether to create starting topography model. If True, must provide starting_topography_kwargs, if False must provide starting_topography by default False
create_starting_prisms (bool, optional) – Choose whether to create starting prisms model. If False, must provide prisms model, by default False
calculate_starting_gravity (bool, optional) – Choose whether to calculate starting gravity from prisms model. If False, must provide column “starting_gravity” in grav_df , by default False
calculate_regional_misfit (bool, optional) – Choose whether to calculate regional misfit. If False, must provide column “reg” in grav_df, if True, must provide`regional_grav_kwargs`, by default False
run_damping_cv (bool, optional) – Choose whether to run cross validation for damping, if True, must supplied damping values with kwarg damping_limits, by default False
run_zref_or_density_cv (bool, optional) – Choose whether to run cross validation for zref or density, if True, must provide zref values, density values, or both with kwargs zref_values or ` density_values`, by default False
run_kfolds_zref_or_density_cv (bool, optional) – Choose whether to run internal kfolds cross validation for zref or density, if True, must provide kwargs for splitting constraints into test/train points with split_kwargs, by default False
plot_cv (bool, optional) – Choose whether to plot the cross validation results, by default False
fname (str, optional) – filename and path to use for saving results. If running a damping CV, will save the study to <fname>_damping_cv_study.pickle and the tuple of the best inversion results to <fname>_damping_cv_results.pickle. If running a density/zref CV, will save the study to <fname>_zref_density_cv_study.pickle and the tuple of the best inversion results to <fname>_zref_density_cv_results.pickle. The final inversion result for all methods will be saved to <fname>_results.pickle, by default will be “tmp_<x>” where x is a random integer between 0 and 999.
kwargs (Any) – keyword arguments for the workflow and inversion, such as starting_topography_kwargs, regional_grav_kwargs, split_kwargs and all the other kwargs supplied to run_inversion.
- Returns:
prisms_df (pandas.DataFrame) – prism properties for each iteration, including columns starting_topo for the unchanged starting topography model and topo for the updated topography model
gravity (pandas.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[pandas.DataFrame, pandas.DataFrame, dict[str, Any], float]