invert4geom.inversion ===================== .. py:module:: invert4geom.inversion Functions --------- .. autoapisummary:: invert4geom.inversion.grav_column_der invert4geom.inversion.jacobian_annular invert4geom.inversion._prism_properties invert4geom.inversion.jacobian_prism invert4geom.inversion.jacobian invert4geom.inversion.solver invert4geom.inversion.update_l2_norms invert4geom.inversion.end_inversion invert4geom.inversion.update_gravity_and_misfit invert4geom.inversion.run_inversion invert4geom.inversion.run_inversion_workflow Module Contents --------------- .. py:function:: 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 :footcite:p:`mccubbineairborne2016`. :param grav_easting: coordinates of gravity observation points. :type grav_easting: numpy.ndarray :param grav_northing: coordinates of gravity observation points. :type grav_northing: numpy.ndarray :param grav_upward: coordinates of gravity observation points. :type grav_upward: numpy.ndarray :param prism_easting: coordinates of prism's center in northing, easting, and upward directions, respectively :type prism_easting: numpy.ndarray :param prism_northing: coordinates of prism's center in northing, easting, and upward directions, respectively :type prism_northing: numpy.ndarray :param prism_top: coordinates of prism's center in northing, easting, and upward directions, respectively :type prism_top: numpy.ndarray :param prism_spacing: resolution of prism layer in meters :type prism_spacing: float :param prism_density: density of prisms, in kg/m^3 :type prism_density: numpy.ndarray :returns: array of vertical derivative of gravity at observation point for series of prisms :rtype: numpy.ndarray .. rubric:: References .. footbibliography:: .. py:function:: 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 :footcite:p:`mccubbineairborne2016`, and calculates it's vertical gravity derivative. Takes arrays from `jacobian`, feeds them into `grav_column_der`, and returns the jacobian. :param grav_easting: coordinates of gravity observation points :type grav_easting: numpy.ndarray :param grav_northing: coordinates of gravity observation points :type grav_northing: numpy.ndarray :param grav_upward: coordinates of gravity observation points :type grav_upward: numpy.ndarray :param prism_easting: coordinates of prism's center in northing, easting, and upward directions, respectively :type prism_easting: numpy.ndarray :param prism_northing: coordinates of prism's center in northing, easting, and upward directions, respectively :type prism_northing: numpy.ndarray :param prism_top: coordinates of prism's center in northing, easting, and upward directions, respectively :type prism_top: numpy.ndarray :param prism_density: density of prisms, in kg/m^3 :type prism_density: numpy.ndarray :param prism_spacing: resolution of prism layer in meters :type prism_spacing: float :param jac: empty jacobian matrix with a row per gravity observation and a column per prism :type jac: numpy.ndarray :returns: returns a jacobian matrix of shape (number of gravity points, number of prisms) :rtype: numpy.ndarray .. rubric:: References .. footbibliography:: .. py:function:: _prism_properties(prisms_layer, method = 'itertools') extract prism properties from prism layer :param prisms_layer: harmonica prism layer :type prisms_layer: xarray.Dataset :param method: choice of method to extract properties, by default "itertools" :type method: str, optional :returns: array of prism properties :rtype: numpy.ndarray .. py:function:: 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. :param prisms_properties: 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 :type prisms_properties: numpy.ndarray :param grav_easting: coordinates of gravity observation points. :type grav_easting: numpy.ndarray :param grav_northing: coordinates of gravity observation points. :type grav_northing: numpy.ndarray :param grav_upward: coordinates of gravity observation points. :type grav_upward: numpy.ndarray :param delta: thickness in meters of small prisms used to calculate vertical derivative :type delta: float :param jac: empty jacobian matrix with a row per gravity observation and a column per prism :type jac: numpy.ndarray :returns: returns a numpy.ndarray of shape (number of gravity points, number of prisms) :rtype: numpy.ndarray .. py:function:: 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 :param deriv_type: choose between "annulus" and "prisms" methods of calculating the vertical derivative of gravity of a prism :type deriv_type: str :param coordinates: coordinate dataframe of gravity observation points with columns "easting", "northing", "upward" :type coordinates: pandas.DataFrame :param empty_jac: optionally provide an empty jacobian matrix of shape (number of gravity observations x number of prisms), by default None :type empty_jac: numpy.ndarray, optional :param prisms_layer: harmonica prism layer, by default None :type prisms_layer: xarray.Dataset, optional :param prism_spacing: resolution of prism layer, by default None :type prism_spacing: float, optional :param prism_size: height of prisms for small prism vertical derivative method, by default None :type prism_size: float, optional :param prisms_properties_method: method for extracting prism properties, by default "itertools" :type prisms_properties_method: str, optional :returns: a filled out jacobian matrix :rtype: numpy.ndarray .. py:function:: 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 :param jac: input jacobian matrix with a row per gravity observation, and a column per prisms. :type jac: numpy.ndarray :param residuals: array of gravity residuals :type residuals: numpy.ndarray :param damping: positive damping (Tikhonov 0th order) regularization :type damping: float | None, optional :param solver_type: choose which solving method to use, by default "scipy least squares" :type solver_type: str, optional :returns: array of correction values to apply to each prism. :rtype: numpy.ndarray .. py:function:: update_l2_norms(current_rmse, last_l2_norm) 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%. :param current_rmse: root mean square error of the residual gravity misfit :type current_rmse: float :param last_l2_norm: l2 norm of the residual gravity misfit :type last_l2_norm: float :returns: * **updated_l2_norm** (*float*) -- the updated l2 norm * **updated_delta_l2_norm** (*float*) -- the updated delta l2 norm .. py:function:: 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) check if the inversion should be terminated :param iteration_number: the iteration number, starting at 1 not 0 :type iteration_number: int :param max_iterations: the maximum allowed iterations, inclusive and starting at 1 :type max_iterations: int :param l2_norms: a list of each iteration's l2 norm :type l2_norms: float :param l2_norm_tolerance: the l2 norm value to end the inversion at :type l2_norm_tolerance: float :param delta_l2_norm: the current iteration's delta l2 norm :type delta_l2_norm: float :param previous_delta_l2_norm: the delta l2 norm of the previous iteration :type previous_delta_l2_norm: float :param delta_l2_norm_tolerance: the delta l2 norm value to end the inversion at :type delta_l2_norm_tolerance: float :param perc_increase_limit: the set tolerance for decimal percentage increase relative to the starting l2 norm :type perc_increase_limit: float :returns: * **end** (*bool*) -- whether or not to end the inversion * **termination_reason** (*list[str]*) -- a list of termination reasons .. py:function:: update_gravity_and_misfit(gravity_df, prisms_ds, 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. :param gravity_df: gravity dataframe with gravity observation coordinate columns ('easting', 'northing'), a gravity data column 'gravity_anomaly', and a regional gravity column ('reg'). :type gravity_df: pandas.DataFrame :param prisms_ds: harmonica prism layer :type prisms_ds: xarray.Dataset :param iteration_number: iteration number to use in updated column names :type iteration_number: int :returns: a gravity dataframe with 2 new columns, one for the iterations forward gravity and one for the iterations residual misfit. :rtype: pandas.DataFrame .. py:function:: 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) 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` :param grav_df: dataframe with gravity data and coordinates, must have columns "gravity_anomaly", "misfit", "res" and "reg", and coordinate columns "easting" and "northing". :type grav_df: pandas.DataFrame :param prism_layer: starting prism layer :type prism_layer: xarray.Dataset :param max_iterations: the maximum allowed iterations, inclusive and starting at 1 :type max_iterations: int :param l2_norm_tolerance: _the l2 norm value to end the inversion at, by default 0.2 :type l2_norm_tolerance: float, optional :param delta_l2_norm_tolerance: the delta l2 norm value to end the inversion at, by default 1.001 :type delta_l2_norm_tolerance: float, optional :param perc_increase_limit: the set tolerance for decimal percentage increase relative to the starting l2 norm, by default 0.10 :type perc_increase_limit: float, optional :param deriv_type: either "annulus" or "prism" to determine method of calculating the vertical derivate of gravity of a prism, by default "annulus" :type deriv_type: str, optional :param jacobian_prism_size: height of prisms in meters for vertical derivative, by default 1 :type jacobian_prism_size: float, optional :param solver_type: solver type to use, by default "scipy least squares" :type solver_type: str, optional :param solver_damping: damping parameter for regularization of the solver, by default None :type solver_damping: float | None, optional :param upper_confining_layer: topographic layer to use as upper limit for inverted topography, by default None :type upper_confining_layer: xarray.DataArray | None, optional :param lower_confining_layer: topographic layer to use as lower limit for inverted topography, by default None :type lower_confining_layer: xarray.DataArray | None, optional :param apply_weighting_grid: use "weighting_grid" to scale surface corrections grid, by default False, by default False :type apply_weighting_grid: bool, optional :param plot_convergence: plot the misfit convergence, by default False :type plot_convergence: bool, optional :param plot_dynamic_convergence: plot the misfit convergence dynamically, by default False :type plot_dynamic_convergence: bool, optional :param results_fname: filename to save results to, by default None :type results_fname: str, optional :param progressbar: set to False to turn off the inversion progressbar :type progressbar: bool, optional :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 .. py:function:: 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) 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 :param grav_df: 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'. :type grav_df: pandas.DataFrame :param create_starting_topography: Choose whether to create starting topography model. If True, must provide `starting_topography_kwargs`, if False must provide `starting_topography` by default False :type create_starting_topography: bool, optional :param create_starting_prisms: Choose whether to create starting prisms model. If False, must provide prisms model, by default False :type create_starting_prisms: bool, optional :param calculate_starting_gravity: Choose whether to calculate starting gravity from prisms model. If False, must provide column "starting_gravity" in grav_df , by default False :type calculate_starting_gravity: bool, optional :param calculate_regional_misfit: 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 :type calculate_regional_misfit: bool, optional :param run_damping_cv: Choose whether to run cross validation for damping, if True, must supplied damping values with kwarg `damping_limits`, by default False :type run_damping_cv: bool, optional :param run_zref_or_density_cv: 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 :type run_zref_or_density_cv: bool, optional :param run_kfolds_zref_or_density_cv: 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 :type run_kfolds_zref_or_density_cv: bool, optional :param plot_cv: Choose whether to plot the cross validation results, by default False :type plot_cv: bool, optional :param fname: filename and path to use for saving results. If running a damping CV, will save the study to _damping_cv_study.pickle and the tuple of the best inversion results to _damping_cv_results.pickle. If running a density/zref CV, will save the study to _zref_density_cv_study.pickle and the tuple of the best inversion results to _zref_density_cv_results.pickle. The final inversion result for all methods will be saved to _results.pickle, by default will be "tmp_" where x is a random integer between 0 and 999. :type fname: str, optional :param kwargs: 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`. :type kwargs: typing.Any :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