Code Monkey home page Code Monkey logo

harmonica's Introduction

Fatiando a Terra


Website | Documentation | Gallery | Docs (development version) | Mailing list

Latest version on PyPI

Travis CI build status

AppVeyor build status

Test coverage status

Code health report by landscape.io

doi:10.5281/zenodo.157746

Disclaimer

Fatiando is under active development and we are still changing the API between releases. Names will change and functions will move as we improve our design. You might have to update your scripts and notebooks to get the latest features from a new release.

Please bear with us.

Overview

Our goal is provide a comprehensive and extensible framework for geophysical data analysis and the development of new methodologies.

Research: Make your research more reproducible by writing a Python script or Jupyter notebook instead of clicking through complicated menus.

Development: Don't start from scratch! Build upon the existing tools in Fatiando to develop new methods.

Teaching: Combine Fatiando with the Jupyter notebook to make rich, interactive documents. Great for teaching fundamental concepts of geophysics!

Getting started

  1. Install Fatiando and its dependencies.
  2. Browse the Gallery for examples of what Fatiando can do.
  3. Take a look at the rest of the Documentation for more information about the library.
  4. Get involved in the community and see how you can help in our Contributor Guide.

Contributing and asking for help

Subscribe to our Google Groups mailing list to stay informed and ask for help: groups.google.com/d/forum/fatiando

We'll post updates to the list about new releases and features, events, and future plans for the project. Get involved to help us shape the project and make it even better!

Another option for reaching out and reporting bugs is to open an issue on Github.

We have an open development process where everything is discussed through Github issues. Anyone can comment and give feedback. See our Roadmap for v1.0 to get a feeling for where the project is headed. Your input is welcome!

Supporting

If you use Fatiando in your research, please cite it in your publications as:

Uieda, L., V. C. Oliveira Jr, and V. C. F. Barbosa (2013), Modeling the Earth with Fatiando a Terra, Proceedings of the 12th Python in Science Conference, pp. 91 - 98.

Please also cite the method papers of individual functions/classes. References are available in the documentation of each module/function/class.

See the CITATION.rst file or the Citing section of the docs for more information.

You can also show your support by buying a sticker from Stickermule. We don't make any money from the sales but it helps spread the word about the project.

License

Fatiando a Terra is free software: you can redistribute it and/or modify it under the terms of the BSD 3-clause License. A copy of this license is provided in LICENSE.txt.

harmonica's People

Contributors

aguspesce avatar benjmy avatar birocoles avatar dependabot[bot] avatar esteban82 avatar fatiando-bot avatar indiauppal avatar leomiquelutti avatar leouieda avatar ll-geo avatar mdtanker avatar mgomezn avatar nshea3 avatar richardscottoz avatar rowanc1 avatar santisoler avatar yagomcastro avatar ziebam avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

harmonica's Issues

Normal gravity calculation

Description of the desired feature

The normal gravity calculation should be done using the Li and Goetze (2001) formula that allows computation at any height. It would use the current ellipsoid (harmonica.get_ellipsoid) to do the calculations. This is what the function signature should look like:

def normal_gravity(latitude, height):
    ...
    return result

Latitude and heights can be any array like object so it's best to do the calculations with numpy functions instead of math.

It might be a bit tricky to test. One way is to compare the computed values at the equator and pole with the ones provided by the ReferenceEllipsoid (which are tested against Hofmann-Wellenhof & Moritz, 2006). That should probably be enough.

Decorator require_numba does not work as expected

Description of the problem

The require_numba decorator for test functions does not work as expected. It decorates the function but the test function itself is not run, so pytest assumes it passed although it has never been run.

This is due to the fact that the function wrapper defined inside require_numba does not run the function passed as argument.

Add sample station gravity data

Description of the desired feature

We have a global gravity grid from spherical harmonics, which is great to show gravity corrections and global methods. It would be great to have some more local station data to show interpolations, localized inversion, etc.

There are a lot of datasets at https://www.ngdc.noaa.gov/mgg/gravity/

One that I find promising is this South African collection: https://www.ngdc.noaa.gov/mgg/gravity/1999/data/regional/africa/

It has different spatial coverages, which is great for showing and testing interpolation, and covers the Bushveld, which might have some interesting anomalies for modeling.

We should add this to our sample data. Probably will need to have a look first and see if there is anything in there or what the data quality is like. We might have to use a subset of the data if the files get too large. It will also probably require a geoid model since heights are relative to sea level. It would probably be good enough to interpolate from a global grid (#29).

Our topography data should probably be referenced to the geoid

Description of the problem

Currently, our sample topography data based on ETOPO1 was converted to be referenced to the WGS84 ellipsoid by adding a geoid model based on EIGEN-6c4. For some applications, like Bouguer corrections on land, this is best. But Bouguer corrections in the oceans require knowing the geoid to account for the water mass between the ellipsoid and the geoid. This also plays a role in Airy isostatic calculations (so it might be best to use geoid reference topography for this).

Since ETOPO1 is originally referenced to "mean sea level", we should probably keep it as is and include a geoid grid as well. Docstrings for the fetch_* functions should be explicit as to which reference is used. Examples of calculating Bouguer anomalies or Airy isostatic roots should also be explicit in which heights they are using and why.

Gallery example of gravity disturbance calculation (local)

We have an example of global gravity disturbance calculation using spherical harmonic data. This is a global grid and is straight forward to calculate the geometric heights using gridded geoid values. It would be great to have a similar example using our land gravity data (#99). For this, we would have to interpolate the geoid values on the observation points (using verde.ScipyGridder with linear interpolation should be enough) so we can get geometric (ellipsoid) heights.

Implement the fast iterative EQL interpolator

Description of the desired feature

Would be nice to implement the fast iterative equivalent layer technique described on Siqueira, Oliveira Jr. and Barbosa (2017). It could be a fast alternative to the EQLHarmonic gridder from #78 that could be very useful when working with a very high number of data points.

The new class should inherit from verde.base.BaseGridder and follow its protocols. See the EQLHarmonic class defined on #78 to get inspiration.

Are you willing to help implement and maintain this feature? Yes

Point masses gravity computation on Cartesian coordinates

Description of the desired feature

The most simple geometry that generates gravitational fields are point masses.
We should make a similar implementation that the one proposed on #47 but for Cartesian coordinates.

Are you willing to help implement and maintain this feature? Yes/No

Move the ReferenceEllipsoid code to Boule

In discussions with @MarkWieczorek from SHTools, we found that there is a lot of repeated code between the projects for handling ellipsoids and their parameters. There are also other packages there that might benefit from this code. So it makes sense to split this out of Harmonica and into a dedicated package which we can all use and develop to keep the dependency light.

I created Boule to serve as this package. The repository has all of the infrastructure used by Fatiando projects and we should move the relevant ellipsoid code there (including the coordinate conversions). Normal gravity calculations should probably remain here, though. But I'm open to discussions on this.

This way, Harmonica remains a strickly gravity and magnetics package.

The main challenge will be retaining the global ellipsoid setting functionality. It will probably need to be a bit more sophisticated and allow setting package-specific ellipsoids. For example:

import boule
import harmonica as hm

boule.set_ellipsoid(boule.GRS80, packages=["harmonica"])

gamma = hm.normal_gravity(...)

Otherwise, setting an ellipsoid for Harmonica would influence SHTools indirectly. Which we might not want since it adds difficult to test side effects. We could have a "global" setting but I usually prefer to set things explicitly so everyone knows what is happening by reading the code.

Enhance EQL gridder to interpolate data in spherical coordinates

Description of the desired feature

Would be nice to improve the EQL Harmonic gridder to be able to interpolate data in spherical coordinates (without the need to project it). This feature would be very useful when gridding potential data on regional or continental scales where the curvature of the Earth must be taken into account. It can be done by modifying the EQLHarmonic class, allowing to pass a coordinate_system argument which will regulate how to compute the Green's function (1 / distance) between the point sources and the observation points.

Are you willing to help implement and maintain this feature? Yes

Standarize Cartesian coordinates with z pointing up

Description of the desired feature

After a discussion, @leouieda and I agreed that would be better to standarize the vertical coordinate in Cartesian coordinate system pointing upwards.
The benefits of this decision are:

  • Computational heights would be positive values, saving users from converting their height coordinates of observation points to negative values.
  • The same goes for structures above the zero height, like topography. Having negative topography heights is not intuitive.
  • When specifying some source geometries (like prisms), their boundaries can be standardized as the ones in spherical coordinates are defined: w, e, s, n, bottom, top. If the vertical direction points upwards, the bottom value will be numerically lower than the top one, in the same way as tesseroids are defined.

One of the drawbacks of this idea is that the coordinate system would be left-handed. Nevertheless, this can be solved with good documentation and careful implementation of the math beneath the forward models.

This new idea also change how gravity acceleration is defined. For sake of standarization, we propose that the forward models that computes the vertical (or radial) component of the gravity acceleration return the downward component rather than the one oriented according to the new vertical direction. We agreed on this because among geophysicists the gravity is always expected to be pointing downward, thus a positive density contrast produces a positive anomaly.
This should be well specified on the documentation and on the future page about coordinates systems (#68).

Are you willing to help implement and maintain this feature? Yes

Raise error if masses array has not the same number of elements as points

Description of the desired feature

On the point mass forward model would be nice to check if the masses has the same number of elements as the points coordinates inside points. On the tesseroid forward model this is done, but for some reason we forgot to add it on the point mass case.

Are you willing to help implement and maintain this feature? Yes

Tesseroid forward modeling

Description of the desired feature

We need to port the tesseroid modeling method (with constant density) to Harmonica. I'm prone to trying numba again because it's a lot more stable now than it was a few years ago. Plus, this will allow us to expand to variable density tesseroids with more ease than Cython. The user can provide a JIT compiled density function to get fast results, which can't be done easily in Cython.

Airy Isostasy calculations

Description of the desired feature

We need a function that calculates the Moho depth from topography based on an Airy isostatic model. This will be useful for calculating isostatic gravity anomalies, for example.

The function only needs to know the topographic height values and crust, mantle, and water densities (optional):

def isostasy_airy(topography, density_crust, density_mantle, density_water=None):
    ...
    return moho_depth

If water density is not given, negative topography will be ignored and the Moho depth for those values will be set to NaN.

The function can live in a new isostasy.py module.

I'm entirely happy with the function name. Any suggestions?

Add function to create tesseroid layers

Description of the desired feature

A tesseroid layer could be defined as a collection of tesseroids sorted in a regular way within a region, whose longitudinal and latitudinal boundaries touch each other leaving no empty space in between nor overlapping each other.
Tesseroids layers are very useful for gravity inversion (e.g. Moho inversion). So the function should allow us to modify the top and bottom boundaries of the tesseroid layer according to a given topography data or to a predicted radii by an inversion algorithm.

The tesseroid layer should be returned as an xarray object, what will ease the process of computing the forward model of the layer.

Are you willing to help implement and maintain this feature? Yes

Repeated points on sample datasets downloaded from ICGEM

Description of the problem

We currently have three sample datasets that have been downloaded from the ICGEM Calculation Service:

  • Earth Topography
  • Earth Gravity
  • Earth Geoid

Each of these models has been created on a regular grid in a geodetic coordinate system.
On the latitude direction they span between -90 and 90 degrees, while on the longitudinal one they do so between -180 and 180 degrees.
The problem is that the -180 and 180 longitudinal degrees are equivalent, therefore each dataset contains one row of repeated values.

Would be nice to regenerate and upload these datasets without any repeated values.
What has to be done?

  1. Fetch each dataset with the existing functions (e.g. hm.datasets.fetch_topography_earth()).
  2. Remove the repeated longitudinal coordinate and its corresponding data values.
  3. Save the new dataset to a NetCDF file (.nc)
  4. Compress it to a .xz file.
  5. Replace the existing files on data with the new ones.
  6. Get the sha256 hash of the new file and update harmonica/datasets/registry.txt

Register pytest use_numba mark

Description of the desired feature

When creating pytest custom marks, they must be registered. Otherwise a PytestUnknownMarkWarning is raised in order to catch possible typos.
We must register use_numba mark according to instructions in https://docs.pytest.org/en/latest/mark.html.

Are you willing to help implement and maintain this feature? Yes

Add variable density tesseroids forward modelling

Description of the desired feature

Would be nice to implement the forward modelling method to compute the gravitational fields of tesseroids with variable densities described on Soler et al. (2019).

The main thing to do is adding the distance-based discretization algorithm, which should be implemented inside Numba jitted functions without allocating memory.
Also, the tesseroid_gravity function should take a density function per tesseroid as argument, which should be jitted in order to speed things up.

Are you willing to help implement and maintain this feature? Yes

Remove AppVeyor CI

Description of the desired feature

It's very slow and we don't need it anymore for testing on Windows now that we have Azure working. To disable it, we just need to remove the .appveyor.yml file from the repository.

Interger values in sample grids should be converted to float

Description of the desired feature

The height in our sample gravity and topography are stored as integers to save space. However, this can cause problems in computations involving division. To avoid having to explain this to users or add warnings, we should convert these data types to float32 in the fetch_* functions once they are loaded.

Point masses gravity computation on geographic coordinates

Description of the desired feature

The most simple geometry that generates gravitational fields are point masses.
Specially in spherical coordinates (or geodetic), where they represent one of the easier to implement forward models.
They can also be used to grid data through Equivalent Sources algorithms and to build tests for future functions.

Moreover, on #46 we propose to split the Tesseroid forward modelling into the discretization and the forward modelling itself through a GLQ. The latter can be addressed by the computation of the gravitational fields generated by point masses. So the point masses implementation would serve as basis for the Tesseroid forward modelling.

Are you willing to help implement and maintain this feature? Yes

Numba deprecation of reflected list

Description of the problem

Pytest is currently raising a NumbaPendingDeprecationWarning due to a scheduled deprecation of type reflected list.
I believe it's raised when a list is passed as argument to a jitted function.

Full code that generated the error

make test

Full error message

# Run a tmp folder to make sure the tests are run on the installed version
mkdir -p tmp-test-dir-with-unique-name
cd tmp-test-dir-with-unique-name; NUMBA_DISABLE_JIT=1 MPLBACKEND='agg' pytest --cov-config=../.coveragerc --cov-report=term-missing --cov=harmonica --doctest-modules -v --pyargs harmonica
====================================== test session starts ======================================
platform linux -- Python 3.7.3, pytest-5.0.1, py-1.8.0, pluggy-0.12.0 -- /home/santi/.miniconda3/envs/harmonica/bin/python
cachedir: .pytest_cache
rootdir: /home/santi/documentos/geo/git/harmonica
plugins: cov-2.7.1
collected 72 items                                                                              

../harmonica/coordinates.py::harmonica.coordinates.geodetic_to_spherical PASSED           [  1%]
../harmonica/coordinates.py::harmonica.coordinates.spherical_to_geodetic PASSED           [  2%]
../harmonica/ellipsoid.py::harmonica.ellipsoid.ReferenceEllipsoid PASSED                  [  4%]
../harmonica/ellipsoid.py::harmonica.ellipsoid.get_ellipsoid PASSED                       [  5%]
../harmonica/ellipsoid.py::harmonica.ellipsoid.print_ellipsoids PASSED                    [  6%]
../harmonica/ellipsoid.py::harmonica.ellipsoid.set_ellipsoid PASSED                       [  8%]
../harmonica/forward/tesseroid.py::harmonica.forward.tesseroid.tesseroid_gravity PASSED   [  9%]
../harmonica/tests/test_coordinates.py::test_geodetic_to_spherical_with_spherical_ellipsoid PASSED [ 11%]
../harmonica/tests/test_coordinates.py::test_geodetic_to_spherical_on_equator PASSED      [ 12%]
../harmonica/tests/test_coordinates.py::test_geodetic_to_spherical_on_poles PASSED        [ 13%]
../harmonica/tests/test_coordinates.py::test_spherical_to_geodetic_with_spherical_ellipsoid PASSED [ 15%]
../harmonica/tests/test_coordinates.py::test_spherical_to_geodetic_on_equator PASSED      [ 16%]
../harmonica/tests/test_coordinates.py::test_spherical_to_geodetic_on_poles PASSED        [ 18%]
../harmonica/tests/test_coordinates.py::test_spherical_to_geodetic_identity PASSED        [ 19%]
../harmonica/tests/test_ellipsoid.py::test_set_ellipsoid PASSED                           [ 20%]
../harmonica/tests/test_ellipsoid.py::test_set_ellipsoid_by_name PASSED                   [ 22%]
../harmonica/tests/test_ellipsoid.py::test_set_ellipsoid_by_name_context PASSED           [ 23%]
../harmonica/tests/test_ellipsoid.py::test_ellipsoid_context PASSED                       [ 25%]
../harmonica/tests/test_gravity_corrections.py::test_normal_gravity PASSED                [ 26%]
../harmonica/tests/test_gravity_corrections.py::test_normal_gravity_arrays PASSED         [ 27%]
../harmonica/tests/test_gravity_corrections.py::test_normal_gravity_non_zero_height PASSED [ 29%]
../harmonica/tests/test_gravity_corrections.py::test_bouguer_correction PASSED            [ 30%]
../harmonica/tests/test_gravity_corrections.py::test_bouguer_correction_zero_topo PASSED  [ 31%]
../harmonica/tests/test_gravity_corrections.py::test_bouguer_correction_xarray PASSED     [ 33%]
../harmonica/tests/test_icgem.py::test_load_icgem_gdf PASSED                              [ 34%]
../harmonica/tests/test_icgem.py::test_load_icgem_gdf_with_height PASSED                  [ 36%]
../harmonica/tests/test_icgem.py::test_load_icgem_gdf_usecols PASSED                      [ 37%]
../harmonica/tests/test_icgem.py::test_missing_shape PASSED                               [ 38%]
../harmonica/tests/test_icgem.py::test_missing_size PASSED                                [ 40%]
../harmonica/tests/test_icgem.py::test_corrupt_shape PASSED                               [ 41%]
../harmonica/tests/test_icgem.py::test_missing_cols_names PASSED                          [ 43%]
../harmonica/tests/test_icgem.py::test_missing_units PASSED                               [ 44%]
../harmonica/tests/test_icgem.py::test_missing_empty_line PASSED                          [ 45%]
../harmonica/tests/test_icgem.py::test_missing_attribute PASSED                           [ 47%]
../harmonica/tests/test_icgem.py::test_missing_lat_lon_attributes PASSED                  [ 48%]
../harmonica/tests/test_icgem.py::test_diff_attrs_vs_cols PASSED                          [ 50%]
../harmonica/tests/test_icgem.py::test_missing_area PASSED                                [ 51%]
../harmonica/tests/test_icgem.py::test_corrupt_area PASSED                                [ 52%]
../harmonica/tests/test_icgem.py::test_empty_file PASSED                                  [ 54%]
../harmonica/tests/test_isostasy.py::test_isostasy_airy_zero_topography PASSED            [ 55%]
../harmonica/tests/test_isostasy.py::test_isostasy_airy PASSED                            [ 56%]
../harmonica/tests/test_isostasy.py::test_isostasy_airy_dataarray PASSED                  [ 58%]
../harmonica/tests/test_point_mass.py::test_invalid_field PASSED                          [ 59%]
../harmonica/tests/test_point_mass.py::test_point_mass_on_origin PASSED                   [ 61%]
../harmonica/tests/test_point_mass.py::test_point_mass_same_radial_direction PASSED       [ 62%]
../harmonica/tests/test_point_mass.py::test_point_mass_potential_on_equator PASSED        [ 63%]
../harmonica/tests/test_point_mass.py::test_point_mass_potential_on_same_meridian PASSED  [ 65%]
../harmonica/tests/test_sample_data.py::test_geoid_earth PASSED                           [ 66%]
../harmonica/tests/test_sample_data.py::test_gravity_earth PASSED                         [ 68%]
../harmonica/tests/test_sample_data.py::test_topography_earth PASSED                      [ 69%]
../harmonica/tests/test_sample_data.py::test_rio_magnetic PASSED                          [ 70%]
../harmonica/tests/test_tesseroid.py::test_single_tesseroid PASSED                        [ 72%]
../harmonica/tests/test_tesseroid.py::test_invalid_field PASSED                           [ 73%]
../harmonica/tests/test_tesseroid.py::test_invalid_distance_size_ratii PASSED             [ 75%]
../harmonica/tests/test_tesseroid.py::test_invalid_density_array PASSED                   [ 76%]
../harmonica/tests/test_tesseroid.py::test_invalid_tesseroid PASSED                       [ 77%]
../harmonica/tests/test_tesseroid.py::test_point_inside_tesseroid PASSED                  [ 79%]
../harmonica/tests/test_tesseroid.py::test_stack_overflow_on_adaptive_discretization PASSED [ 80%]
../harmonica/tests/test_tesseroid.py::test_distance_tesseroid_point PASSED                [ 81%]
../harmonica/tests/test_tesseroid.py::test_tesseroid_dimensions PASSED                    [ 83%]
../harmonica/tests/test_tesseroid.py::test_longitude_continuity PASSED                    [ 84%]
../harmonica/tests/test_tesseroid.py::test_longitude_continuity_equivalent_tesseroids PASSED [ 86%]
../harmonica/tests/test_tesseroid.py::test_split_tesseroid PASSED                         [ 87%]
../harmonica/tests/test_tesseroid.py::test_split_tesseroid_only_longitude PASSED          [ 88%]
../harmonica/tests/test_tesseroid.py::test_split_tesseroid_only_latitude PASSED           [ 90%]
../harmonica/tests/test_tesseroid.py::test_split_tesseroid_only_radius PASSED             [ 91%]
../harmonica/tests/test_tesseroid.py::test_split_tesseroid_only_horizontal PASSED         [ 93%]
../harmonica/tests/test_tesseroid.py::test_adaptive_discretization_on_radii SKIPPED       [ 94%]
../harmonica/tests/test_tesseroid.py::test_adaptive_discretization_vs_distance_size_ratio SKIPPED [ 95%]
../harmonica/tests/test_tesseroid.py::test_two_dimensional_adaptive_discretization SKIPPED [ 97%]
../harmonica/tests/test_tesseroid.py::test_spherical_shell_two_dim_adaptive_discret SKIPPED [ 98%]
../harmonica/tests/test_tesseroid.py::test_spherical_shell_three_dim_adaptive_discret SKIPPED [100%]

======================================= warnings summary ========================================
/home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/_pytest/mark/structures.py:332
/home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/_pytest/mark/structures.py:332
  /home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/_pytest/mark/structures.py:332: PytestUnknownMarkWarning: Unknown pytest.mark.use_numba - is this a typo?  You can register custom marks to avoid this warning - for details, see https://docs.pytest.org/en/latest/mark.html
    PytestUnknownMarkWarning,

-- Docs: https://docs.pytest.org/en/latest/warnings.html

----------- coverage: platform linux, python 3.7.3-final-0 -----------
Name                                                                         Stmts   Miss Branch BrPart  Cover   Missing
------------------------------------------------------------------------------------------------------------------------
/home/santi/documentos/geo/git/harmonica/harmonica/constants.py                  1      0      0      0   100%
/home/santi/documentos/geo/git/harmonica/harmonica/coordinates.py               32      0      0      0   100%
/home/santi/documentos/geo/git/harmonica/harmonica/datasets/sample_data.py      34      0      0      0   100%
/home/santi/documentos/geo/git/harmonica/harmonica/ellipsoid.py                 68      0     12      0   100%
/home/santi/documentos/geo/git/harmonica/harmonica/forward/point_mass.py        40      0     12      0   100%
/home/santi/documentos/geo/git/harmonica/harmonica/forward/tesseroid.py        177      0     56      0   100%
/home/santi/documentos/geo/git/harmonica/harmonica/gravity_corrections.py       32      0      0      0   100%
/home/santi/documentos/geo/git/harmonica/harmonica/io.py                        65      0     40      0   100%
/home/santi/documentos/geo/git/harmonica/harmonica/isostasy.py                  17      0      2      0   100%
/home/santi/documentos/geo/git/harmonica/harmonica/version.py                    4      0      0      0   100%
------------------------------------------------------------------------------------------------------------------------
TOTAL                                                                          470      0    122      0   100%

======================= 67 passed, 5 skipped, 2 warnings in 8.22 seconds ========================
cp tmp-test-dir-with-unique-name/.coverage* .
rm -rvf tmp-test-dir-with-unique-name
'tmp-test-dir-with-unique-name/.coverage' borrado
'tmp-test-dir-with-unique-name/.coverage.eevaa.16147.400939' borrado
removed directory 'tmp-test-dir-with-unique-name'
# Run a tmp folder to make sure the tests are run on the installed version
mkdir -p tmp-test-dir-with-unique-name
cd tmp-test-dir-with-unique-name; NUMBA_DISABLE_JIT=0 MPLBACKEND='agg' pytest --doctest-modules -v --pyargs -m use_numba harmonica
====================================== test session starts ======================================
platform linux -- Python 3.7.3, pytest-5.0.1, py-1.8.0, pluggy-0.12.0 -- /home/santi/.miniconda3/envs/harmonica/bin/python
cachedir: .pytest_cache
rootdir: /home/santi/documentos/geo/git/harmonica
plugins: cov-2.7.1
collected 72 items / 50 deselected / 22 selected                                                

../harmonica/tests/test_point_mass.py::test_point_mass_on_origin PASSED                   [  4%]
../harmonica/tests/test_point_mass.py::test_point_mass_same_radial_direction PASSED       [  9%]
../harmonica/tests/test_point_mass.py::test_point_mass_potential_on_equator PASSED        [ 13%]
../harmonica/tests/test_point_mass.py::test_point_mass_potential_on_same_meridian PASSED  [ 18%]
../harmonica/tests/test_tesseroid.py::test_single_tesseroid PASSED                        [ 22%]
../harmonica/tests/test_tesseroid.py::test_invalid_tesseroid PASSED                       [ 27%]
../harmonica/tests/test_tesseroid.py::test_point_inside_tesseroid PASSED                  [ 31%]
../harmonica/tests/test_tesseroid.py::test_stack_overflow_on_adaptive_discretization PASSED [ 36%]
../harmonica/tests/test_tesseroid.py::test_distance_tesseroid_point PASSED                [ 40%]
../harmonica/tests/test_tesseroid.py::test_tesseroid_dimensions PASSED                    [ 45%]
../harmonica/tests/test_tesseroid.py::test_longitude_continuity PASSED                    [ 50%]
../harmonica/tests/test_tesseroid.py::test_longitude_continuity_equivalent_tesseroids PASSED [ 54%]
../harmonica/tests/test_tesseroid.py::test_split_tesseroid PASSED                         [ 59%]
../harmonica/tests/test_tesseroid.py::test_split_tesseroid_only_longitude PASSED          [ 63%]
../harmonica/tests/test_tesseroid.py::test_split_tesseroid_only_latitude PASSED           [ 68%]
../harmonica/tests/test_tesseroid.py::test_split_tesseroid_only_radius PASSED             [ 72%]
../harmonica/tests/test_tesseroid.py::test_split_tesseroid_only_horizontal PASSED         [ 77%]
../harmonica/tests/test_tesseroid.py::test_adaptive_discretization_on_radii PASSED        [ 81%]
../harmonica/tests/test_tesseroid.py::test_adaptive_discretization_vs_distance_size_ratio PASSED [ 86%]
../harmonica/tests/test_tesseroid.py::test_two_dimensional_adaptive_discretization PASSED [ 90%]
../harmonica/tests/test_tesseroid.py::test_spherical_shell_two_dim_adaptive_discret PASSED [ 95%]
../harmonica/tests/test_tesseroid.py::test_spherical_shell_three_dim_adaptive_discret PASSED [100%]

======================================= warnings summary ========================================
/home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/_pytest/mark/structures.py:332
/home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/_pytest/mark/structures.py:332
  /home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/_pytest/mark/structures.py:332: PytestUnknownMarkWarning: Unknown pytest.mark.use_numba - is this a typo?  You can register custom marks to avoid this warning - for details, see https://docs.pytest.org/en/latest/mark.html
    PytestUnknownMarkWarning,

harmonica/tests/test_tesseroid.py::test_single_tesseroid
  /home/santi/documentos/geo/git/harmonica/harmonica/forward/tesseroid.py:258: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'glq_nodes' of function 'tesseroids_to_point_masses'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 275:
  @jit(nopython=True)
  def tesseroids_to_point_masses(
  ^
  
    weights,

harmonica/tests/test_tesseroid.py::test_single_tesseroid
  /home/santi/documentos/geo/git/harmonica/harmonica/forward/tesseroid.py:258: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'glq_weights' of function 'tesseroids_to_point_masses'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 275:
  @jit(nopython=True)
  def tesseroids_to_point_masses(
  ^
  
    weights,

harmonica/tests/test_tesseroid.py::test_single_tesseroid
harmonica/tests/test_tesseroid.py::test_single_tesseroid
  /home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/numba/ir_utils.py:1958: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'glq_nodes' of function 'jit_tesseroid_gravity'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 172:
  @jit(nopython=True, parallel=True)
  def jit_tesseroid_gravity(
  ^
  
    warnings.warn(NumbaPendingDeprecationWarning(msg, loc=loc))

harmonica/tests/test_tesseroid.py::test_single_tesseroid
harmonica/tests/test_tesseroid.py::test_single_tesseroid
  /home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/numba/ir_utils.py:1958: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'glq_weights' of function 'jit_tesseroid_gravity'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 172:
  @jit(nopython=True, parallel=True)
  def jit_tesseroid_gravity(
  ^
  
    warnings.warn(NumbaPendingDeprecationWarning(msg, loc=loc))

harmonica/tests/test_tesseroid.py::test_stack_overflow_on_adaptive_discretization
  /home/santi/documentos/geo/git/harmonica/harmonica/forward/tesseroid.py:436: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'coordinates' of function '_distance_tesseroid_point'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 508:
  @jit(nopython=True)
  def _distance_tesseroid_point(
  ^
  
    distance = _distance_tesseroid_point(coordinates, tesseroid)

harmonica/tests/test_tesseroid.py::test_stack_overflow_on_adaptive_discretization
  /home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/numba/ir_utils.py:1958: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'coordinates' of function '_adaptive_discretization'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 380:
  @jit(nopython=True)
  def _adaptive_discretization(
  ^
  
    warnings.warn(NumbaPendingDeprecationWarning(msg, loc=loc))

harmonica/tests/test_tesseroid.py::test_distance_tesseroid_point
  /home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/numba/ir_utils.py:1958: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'coordinates' of function '_distance_tesseroid_point'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 508:
  @jit(nopython=True)
  def _distance_tesseroid_point(
  ^
  
    warnings.warn(NumbaPendingDeprecationWarning(msg, loc=loc))

harmonica/tests/test_tesseroid.py::test_distance_tesseroid_point
  /home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/numba/ir_utils.py:1958: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'tesseroid' of function '_distance_tesseroid_point'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 508:
  @jit(nopython=True)
  def _distance_tesseroid_point(
  ^
  
    warnings.warn(NumbaPendingDeprecationWarning(msg, loc=loc))

harmonica/tests/test_tesseroid.py::test_tesseroid_dimensions
  /home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/numba/ir_utils.py:1958: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'tesseroid' of function '_tesseroid_dimensions'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 492:
  @jit(nopython=True)
  def _tesseroid_dimensions(tesseroid):
  ^
  
    warnings.warn(NumbaPendingDeprecationWarning(msg, loc=loc))

harmonica/tests/test_tesseroid.py::test_longitude_continuity
  /home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/numba/ir_utils.py:1958: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'region' of function '_longitude_continuity'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 574:
  @jit(nopython=True)
  def _longitude_continuity(region):
  ^
  
    warnings.warn(NumbaPendingDeprecationWarning(msg, loc=loc))

harmonica/tests/test_tesseroid.py::test_split_tesseroid
  /home/santi/.miniconda3/envs/harmonica/lib/python3.7/site-packages/numba/ir_utils.py:1958: NumbaPendingDeprecationWarning: 
  Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'tesseroid' of function '_split_tesseroid'.
  
  For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types
  
  File "../harmonica/forward/tesseroid.py", line 467:
  @jit(nopython=True)
  def _split_tesseroid(
  ^
  
    warnings.warn(NumbaPendingDeprecationWarning(msg, loc=loc))

-- Docs: https://docs.pytest.org/en/latest/warnings.html
==================== 22 passed, 50 deselected, 15 warnings in 16.22 seconds =====================
rm -rvf tmp-test-dir-with-unique-name
removed directory 'tmp-test-dir-with-unique-name'

System information

  • Operating system: Manjaro XFCE 64bits
  • Python installation (Anaconda, system, ETS): Miniconda
  • Version of Python: 3.7.3
  • Version of this package: 8222ed7
  • If using conda, paste the output of conda list below:
output of conda list
# packages in environment at /home/santi/.miniconda3/envs/harmonica:
#
# Name                    Version                   Build  Channel
_libgcc_mutex             0.1                        main  
alabaster                 0.7.12                     py_0    conda-forge
appdirs                   1.4.3                      py_1    conda-forge
asn1crypto                0.24.0                py37_1003    conda-forge
astroid                   2.2.5                    py37_0    conda-forge
atomicwrites              1.3.0                      py_0    conda-forge
attrs                     19.1.0                     py_0    conda-forge
babel                     2.7.0                      py_0    conda-forge
black                     19.3b0                     py_0    conda-forge
bzip2                     1.0.6             h14c3975_1002    conda-forge
ca-certificates           2019.6.16            hecc5488_0    conda-forge
cartopy                   0.17.0          py37he1be148_1005    conda-forge
certifi                   2019.6.16                py37_0    conda-forge
cffi                      1.12.3           py37h8022711_0    conda-forge
chardet                   3.0.4                 py37_1003    conda-forge
click                     7.0                        py_0    conda-forge
coverage                  4.5.3            py37h14c3975_0    conda-forge
cryptography              2.7              py37h72c5cf5_0    conda-forge
cycler                    0.10.0                     py_1    conda-forge
dbus                      1.13.6               he372182_0    conda-forge
docutils                  0.14                  py37_1001    conda-forge
entrypoints               0.3                   py37_1000    conda-forge
expat                     2.2.5             he1b5a44_1003    conda-forge
flake8                    3.7.7                    py37_0    conda-forge
fontconfig                2.13.1            he4413a7_1000    conda-forge
freetype                  2.10.0               he983fc9_0    conda-forge
geos                      3.7.2                he1b5a44_1    conda-forge
gettext                   0.19.8.1          hc5be6a0_1002    conda-forge
glib                      2.58.3            h6f030ca_1001    conda-forge
gst-plugins-base          1.14.5               h0935bb2_0    conda-forge
gstreamer                 1.14.5               h36ae1b5_0    conda-forge
harmonica                 0.1.0a1+41.g8222ed7           dev_0    
icu                       58.2              hf484d3e_1000    conda-forge
idna                      2.8                   py37_1000    conda-forge
imagesize                 1.1.0                      py_0    conda-forge
importlib_metadata        0.18                     py37_0    conda-forge
isort                     4.3.21                   py37_0    conda-forge
jinja2                    2.10.1                     py_0    conda-forge
joblib                    0.13.2                     py_0    conda-forge
jpeg                      9c                h14c3975_1001    conda-forge
kiwisolver                1.1.0            py37hc9558a2_0    conda-forge
lazy-object-proxy         1.4.1            py37h516909a_0    conda-forge
libblas                   3.8.0               10_openblas    conda-forge
libcblas                  3.8.0               10_openblas    conda-forge
libffi                    3.2.1             he1b5a44_1006    conda-forge
libgcc-ng                 9.1.0                hdf63c60_0  
libgfortran-ng            7.3.0                hdf63c60_0  
libiconv                  1.15              h516909a_1005    conda-forge
liblapack                 3.8.0               10_openblas    conda-forge
libopenblas               0.3.6                h6e990d7_4    conda-forge
libpng                    1.6.37               hed695b0_0    conda-forge
libstdcxx-ng              9.1.0                hdf63c60_0  
libtiff                   4.0.10            h57b8799_1003    conda-forge
libuuid                   2.32.1            h14c3975_1000    conda-forge
libxcb                    1.13              h14c3975_1002    conda-forge
libxml2                   2.9.9                h13577e0_1    conda-forge
llvmlite                  0.29.0           py37hfd453ef_1    conda-forge
lz4-c                     1.8.3             he1b5a44_1001    conda-forge
markupsafe                1.1.1            py37h14c3975_0    conda-forge
matplotlib                3.1.1                    py37_0    conda-forge
matplotlib-base           3.1.1            py37hfd891ef_0    conda-forge
mccabe                    0.6.1                      py_1    conda-forge
more-itertools            7.1.0                      py_0    conda-forge
ncurses                   6.1               hf484d3e_1002    conda-forge
numba                     0.44.1           py37hb3f55d8_0    conda-forge
numpy                     1.16.4           py37h95a1406_0    conda-forge
numpydoc                  0.9.1                      py_0    conda-forge
olefile                   0.46                       py_0    conda-forge
openblas                  0.3.6                h6e990d7_4    conda-forge
openssl                   1.1.1c               h516909a_0    conda-forge
owslib                    0.18.0                     py_0    conda-forge
packaging                 19.0                       py_0    conda-forge
pandas                    0.24.2           py37hb3f55d8_0    conda-forge
pcre                      8.41              hf484d3e_1003    conda-forge
pillow                    6.1.0            py37he7afcd5_0    conda-forge
pip                       19.1.1                   py37_0    conda-forge
pluggy                    0.12.0                     py_0    conda-forge
pooch                     0.5.2                    py37_0    conda-forge
proj4                     6.1.0                he751ad9_2    conda-forge
pthread-stubs             0.4               h14c3975_1001    conda-forge
py                        1.8.0                      py_0    conda-forge
pycodestyle               2.5.0                      py_0    conda-forge
pycparser                 2.19                     py37_1    conda-forge
pyepsg                    0.4.0                      py_0    conda-forge
pyflakes                  2.1.1                      py_0    conda-forge
pygments                  2.4.2                      py_0    conda-forge
pykdtree                  1.3.1           py37h3010b51_1002    conda-forge
pylint                    2.3.1                    py37_0    conda-forge
pyopenssl                 19.0.0                   py37_0    conda-forge
pyparsing                 2.4.0                      py_0    conda-forge
pyproj                    2.2.1            py37hc44880f_0    conda-forge
pyqt                      5.9.2            py37hcca6a23_0    conda-forge
pyshp                     2.1.0                      py_0    conda-forge
pysocks                   1.7.0                    py37_0    conda-forge
pytest                    5.0.1                    py37_0    conda-forge
pytest-cov                2.7.1                      py_0    conda-forge
python                    3.7.3                h33d41f4_1    conda-forge
python-dateutil           2.8.0                      py_0    conda-forge
pytz                      2019.1                     py_0    conda-forge
qt                        5.9.7                h52cfd70_2    conda-forge
readline                  8.0                  hf8c457e_0    conda-forge
requests                  2.22.0                   py37_0    conda-forge
scikit-learn              0.21.2           py37hcdab131_1    conda-forge
scipy                     1.3.0            py37h921218d_0    conda-forge
setuptools                41.0.1                   py37_0    conda-forge
shapely                   1.6.4           py37hec07ddf_1006    conda-forge
sip                       4.19.8          py37hf484d3e_1000    conda-forge
six                       1.12.0                py37_1000    conda-forge
snowballstemmer           1.9.0                      py_0    conda-forge
sphinx                    1.8.5                    py37_0    conda-forge
sphinx-gallery            0.4.0                    py37_0    conda-forge
sphinx_rtd_theme          0.4.3                      py_0    conda-forge
sphinxcontrib-websupport  1.1.2                      py_0    conda-forge
sqlite                    3.29.0               hcee41ef_0    conda-forge
tk                        8.6.9             hed695b0_1002    conda-forge
toml                      0.10.0                     py_0    conda-forge
tornado                   6.0.3            py37h516909a_0    conda-forge
urllib3                   1.24.3                   py37_0    conda-forge
verde                     1.1.0                      py_0    conda-forge
wcwidth                   0.1.7                      py_1    conda-forge
wheel                     0.33.4                   py37_0    conda-forge
wrapt                     1.11.2           py37h516909a_0    conda-forge
xarray                    0.12.3                     py_0    conda-forge
xorg-libxau               1.0.9                h14c3975_0    conda-forge
xorg-libxdmcp             1.1.3                h516909a_0    conda-forge
xz                        5.2.4             h14c3975_1001    conda-forge
zipp                      0.5.1                      py_0    conda-forge
zlib                      1.2.11            h14c3975_1004    conda-forge
zstd                      1.4.0                h3b9ef0a_0    conda-forge

Move sanity checks of tesseroids outside jitted functions

Description of the desired feature

Would be better to remove the functions that perform sanity checks of tesseroids outside the jitted functions intended to perform the forward modelling. For any arbitrary set of tesseroids, the current location of these functions is not an issue, but for future tesseroid arrangements like layers (#83) or meshes, these functions can be simplified: instead of calling them once for each tesseroid, we could write alternate simpler functions that work on the entire tesseroid collection.

So, the proposal is to refactor:

  • _check_tesseroid(),
  • _check_point_outside_tesseroid() and
  • _longitude_continuity()
    in order to be called from tesseroid_gravity(). We could use some numpy array properties, so there's no need to jit them. They must take the entire set of tesseroids and coordinates as arguments.

Are you willing to help implement and maintain this feature? Yes

Bouguer correction

Description of the desired feature

Our gravity processing will need the Bouguer correction. The correction aims to remove the effect of masses above the reference ellipsoid (topography) or the mass difference below the ellipsoid (oceans). The function should take two densities, one to apply to the masses at height > 0 (above) and another for masses at height < 0 (below). No contrasts should be calculate (i.e., to remove the effect of the oceans the user should give density_below=1040 - 2670). It's OK to provide defaults for above and below densities.

This should be a straight forward function:

def bouguer_plate(height, density_above=2670, density_below=-1630):
    ...
    return results

Notice that this function doesn't perform the correction itself, it calculates the effect of the Bouguer plates that has to be subtracted from the data. This makes it more inline with calculating a topographic effect using prisms or tesseroids.

Add option to avoid sanity checks on forward models

Description of the desired feature

The current forward models for prisms and tesseroids perform sanity checks on the input arrays. Most of them involve looping through computation points and sources, resulting in a task that could take a significant amount of time. Although this time is negligible when performing a single forward computation, it could be problematic when several forward computations are needed, e.g. on inversions. Would be nice to include an optional argument to avoid performing these sanity checks and also add a notice to only disable them when it's really needed and in a very careful way.

A cool name for the new optional argument could be disable_checks that should be False by default.
The sanity checks that should be optionally disabled are the one that make use of a private function, like _check_prisms().

The new optional argument should have a description on the docstring, explaining when you would set it to True and its possible unwanted consequences.

Are you willing to help implement and maintain this feature? Yes

Speed up EQL for high number of data points

Description of the desired feature

Would be nice to improve the harmonic Equivalent Layer gridder in order to speed up its computations when working with a great number of data points.
Firstly, would be nice to make the Jacobian matrix sparse. This could be done by setting to zero those elements that involve distant data and source points.

Further techniques could be related to reducing the computation time involved on the forward modelling, reducing the number of point sources, etc. These alternatives could be investigated separately, so we could focus on the first one.

Are you willing to help implement and maintain this feature? Yes

Missing first eccentricity in ReferenceEllipsoid

Description of the desired feature

The ReferenceEllipsoid class has two eccentricity attributes: linear_eccentricity (distance between the center of the ellipsoid and one of its nodes in meters) and second_eccentricity.
For some calculations, i.e. coordinates conversion, the first eccentricity is needed.
The first and second eccentricity are defined as follows:

first_eccentricity = linear_eccentricty / semimajor_axis
second_eccentricity = linear_eccentricty / semiminor_axis

Although the first eccentricity could be calculated when needed, I think that including it as an attribute improves the user experience, avoiding the missuse of the second eccentricity as it were the first: if the user is not aware of the existence of different eccentricities, he or she may use the second one as it were the first.

Are you willing to help implement and maintain this feature? Yes

Pratt isostasy

Description of the desired feature

PRs #17 and #36 added a function to calculate the Moho depth based on an Airy isostasy model. It would be great to have a Pratt isostasy calculation function as well. It would take the topography/bathymetry, respective densities, and a compensation depth and output the lateral density variation. The function should look something like this:

def isostasy_pratt(topography, density_crust=2800, density_water=1000,
                   compensation_depth=100e3):
    ...
    return density

The body will be very similar to isostasy_airy to allow numpy arrays or xarray.DataArray as input topographies. The docstring would benefit from having a figure to show what the terms in the equations are.

A good example case would be calculate the density variation along the mid-atlantic ridge.

Wrap all docstrings to 79 characters

Description of the desired feature

After fatiando/verde#177 raised by @prisae, fatiando/community#9 and fatiando/community#10 it has been decided that all docstrings must be wrapped to 79 characters per line.
We should warp all the existing docstrings in Harmonica to 79 characters and also configure flake8 to check for this when running the style checks.

Sadly, there's no way (yet) to automatically change all docstrings, but at least we can use flake8 to raise lines that fail. This can be done by setting max-doc-length to 79 characters on the flake8 configuration under setup.cfg.

Replace numpydoc for napoleon to unpin Sphinx

Description of the desired feature

The combination of latest sphinx (2.2.0), numpydoc, and sphinx_rtd_theme
causes a mangled rendering of the parameter list. We were pinning sphinx
to 1.8.5, which is getting old and causing dependency problems.
See discussion in fatiando/pooch#113.

A possible solution is to replace numpydoc with sphinx.ext.napoleon.
It solves the problem although it generates a slightly different output.
The class autosummary template must be modified in order to account for
missing methods summary block.

This replacement has been achieved by @leouieda on fatiando/pooch#122.
See the changes applied on that PR for inspiration.

Function to transform geodetic to geocentric spherical coordinates

Description of the desired feature

The tesseroid modeling is done in geocentric spherical coordinates (longitude, geocentric latitude, radius) but the data is usually in geodetic coordinates (longitude, geodetic latitude, geometric height). Proj4 apparently only converts the latitudes (https://proj4.org/operations/conversions/geoc.html) but not the height to radius. We need a function that does this conversion so that we can use tesseroids appropriately.

The function will use the current ellipsoid (get_ellipsoid) to do the conversion. Tests can use simplified ellipsoids to check that the computations are correct (like setting the flattening to 1 and axis to an integer). It could go into a new harmonica/coordinates.py module and should look like:

def geodetic_to_geocentric(latitude, height):
    ...
    return geocentric_latitude, radius

Unify computation of distances in spherical and Cartesian coordinates

Description of the desired feature

To prevent repeated code, would be nice to have single functions for computing distances between two points in both spherical and Cartesian coordinates. They must be Numba jitted functions to speed things up.

For the spherical coordinates case, the points are given in longitude, latitude and radius coordinates, but several trigonometric functions must be applied in order to ultimately get the distance between these two points. Some forward models precompute these quantities in order to save time. So, would be nice to have a private function that takes all these precomputed quantities as arguments and returns the distance, and a public function that just takes the points' coordinates, compute this quantities and use this private function to get the distance between them.

Moreover, some forward models just needs the distance square, so we could make another private function that just returns the distance square, i.e. doesn't apply the square root.

For example:

@jit(nopython=True)
def distance_spherical(point_a, point_b):
    # Get coordinates of the two points
    longitude, latitude, radius = point_a[:]
    longitude_p, latitude_p, radius_p = point_b[:]
    # Convert angles to radians
    longitude, latitude = np.radians(longitude), np.radians(latitude)
    longitude_p, latitude_p = np.radians(longitude_p), np.radians(latitude_p)
    # Compute trigonometric quantities
    cosphi_p = np.cos(latitude_p)
    sinphi_p = np.sin(latitude_p)
    cosphi = np.cos(latitude)
    sinphi = np.sin(latitude)
    return _distance_spherical(
        longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
    )


@jit(nopython=True)
def _distance_spherical(
    longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
):
    return np.sqrt(
        _distance_sq_spherical(
            longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
        )
    )


@jit(nopython=True)
def _distance_sq_spherical(
    longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
):
    coslambda = np.cos(longitude_p - longitude)
    cospsi = sinphi_p * sinphi + cosphi_p * cosphi * coslambda
    return (radius - radius_p) ** 2 + 2 * radius * radius_p * (1 - cospsi)

Are you willing to help implement and maintain this feature? Yes

Bad code styling when converting radians to degrees

Description of the problem

On geodetic_to_spherical the geocentric_latitude is computed in degrees by converting it by multiplying the factor 180 / numpy.pi. Would be better to use the numpy.degrees() function.

Catch numba errors about mixing int and float in forward modelling

Description of the desired feature

Our forward modelling code (harmonica/forward/) uses numba to speed up computations. It can be a bit picky about input data types and doesn't like it when we mix integers and floats. So specifying a tesseroid with [10, 20, 0, 15, 1000, 1200] will cause an error from numba about mixing int64 and float64. The actual error message is ugly, intimidating, and hard to read for non-experts. It would be better if we caught the exception in a try:... except:... block and added some more user friendly messages to it, explaining that this is probably due to mixing ints and floats in inputs and providing possible solutions (casting input arrays to float, for example).

For testing, it would help to first make a test that passes ints to the functions and verifies that they actually fail (so if numba ever stops raising these errors we can remove the code that catches them).

General equivalent layer solution

Description of the desired feature

The equivalent layer (EQL) methods are great for interpolation, data integration, upward continuation, etc. Our implementations will be based on the Verde gridding classes. The most basic kind would be an EQL that uses 1/r as the basis function. This can be used with any kind of harmonic function but cannot integrate different fields (grav and gradients) and can't reduce to the pole. But it fits the Verde model well. We should implement this as a first step and then make more specialized gridders that can receive multiple fields as input.

The class should inherit from verde.base.BaseGridder and follow its protocols. It could live in a new package called harmonica/equivalent_layer/harmonic.py or something like that.

Sample dataset: Global gravity from ICGEM

Description of the desired feature

A global gravity grid from an ICGEM spherical harmonic model (like EIGEN-6c4) would be great to showcase our gravity processing and any global methods, like tesseroid topographic correction. This should be an xarray.Dataset with raw gravity (no corrections), latitude, longitude, and ellipsoid height. To make things easier, gravity should probably be calculated on top of topography so that Bouguer corrections work reasonably.

Missing attributes on ICGEM datasets

Description of the problem

We currently have three sample datasets that have been downloaded from the ICGEM Calculation Service:

  • Earth Topography
  • Earth Gravity
  • Earth Geoid

Each dataset can be fetched through its own function. Every one of these functions claims to have more details regarding the generation of the grid under the attrs attribute of the Dataset.
The problem is neither of the returned Datasets have non empty attrs.

Would be nice to correct this and include the attributes related to the grid generation.
I think these grids where generated by @leouieda, so maybe he has more idea of what information should be included on the attrs.

Just by naming some possible attributes that appear on the ICGEM Calculation Service:

  • Reference ellipsoid
  • Height above the ellipsoid
  • Grid step
  • Model used (for gravity and geoid)
  • Tide system
  • Zero degree term (yes/no)
  • Gaussian filter and filter length
  • Low-pass filtering
    • Start gentle cut
    • Maximum degree

Add computation of horizontal components of gravitational acceleration

Description of the desired feature

Point mass and tesseroid forward models are currently capable of computing only the potential and the vertical (or radial) component of the gravitational acceleration.
Would be nice to add the ability to compute the horizontal componentes, i.e. the easting and northing ones.

Are you willing to help implement and maintain this feature? Yes

Ugly docstrings generated by Sphinx

Description of the problem

On the dev version of the Harmonica documentation, Sphinx is not compiling the docstrings right.
For example, on normal_gravity docstring the Parameters and Returns mix the variable name with the variable type.

I leave a screenshot of the current normal_gravity docstring.
harmonica-docs

Drop support for Python 3.5

Conda-forge is no longer building 3.5 packages since October 2018. We should probably remove it from our CIs and officially drop support for new releases. This will mean updating:

  • doc/install.rst
  • setup.py
  • .travis.yml and .appveyor.yml

Split up the isostasy code into a separate package

Since I've been rethinking a lot of the Fatiando packages lately, Harmonica is going to escape. If our goal is to be a "gravity and magnetics" package, then isostasy/flexure calculations don't really fit in here. Sure, we use them in conjunction with gravity but actual calculations don't really have any gravity in them.

I propose moving this type of function to a separate isostasy/flexure package. With the new package template, it's a lot easier to create new packages and we can release them as often as we need (without waiting for unrelated things in Harmonica to be resolved). Maintenance is not that complicated anymore and a lot of our process is automated. So maintaining a single large package is not really easier now. This is why this whole split of fatiando is happening in the first place.

As for possible names, the following aren't taken on PyPI:

  • Roots
  • Isostasy
  • Isostatic
  • Jiboia (a play on potuguese words: jiboia is a boa constrictor [snake] and boia means "float" or "bouy")

Avoid casting data type on forward models

Description of the desired feature

The current forward models (prisms, point masses and tesseroids) cast the data type of some (or all) input arrays. The original idea behind it was to avoid working with input arrays full of integers (that could produce wrong results or even Numba errors). Nevertheless, casting the type of the inputs arrays implies copying them, resulting in an unnecessary duplication of memory.
We should stop casting the type of the input arrays on every forward model, and probably raise a warning (or an error) in case the passed arrays are full of integers.

This issue was raised after #100.

Are you willing to help implement and maintain this feature? Yes

Source depth from radial average power spectrum

Description of the desired feature
The Fatiando package provides an example using the Hawaii gravity data and computing the radial average of its power spectrum.
Is it possible to estimate the average source depth according to the radial average of the power spectrum? In fact, there are several papers mentioned this, even the uncertainty is still controversial.
The equation is similar like h=(P(r1)-P(r2))/(4pi(r2-r1))

For a giving 2D anomaly (gravity or magnetic), we could get the radial average of the power spectrum and the possible source depths.

Are you willing to help implement and maintain this feature? Yes/No
Yes

Function to read ICGEM data files

Description of the desired feature

The ICGEM Calculation Service is an amazing resource to get gravity data. Their grids are in an ASCII format that can be tricky to parse because some values are included in the header. It would be great to have a function that reads the ICGEM file and returns an xarray.Dataset.

This will really help with #9

Remove "Forward modelling" from tesseroid gallery example

Description of the desired feature

Because all gallery examples that perform forward modelling now live inside forward, it's redundant to add "Forward modelling" to the example title. Would be nice to remove it and title it as "Tesseroid" or "Single Tesseroid".

Are you willing to help implement and maintain this feature? Yes

Avoid stacking coordinates on tesseroid forward modelling

Description of the desired feature

As raised by @leouieda on #98 and #86, we should avoid stacking coordiantes arrays in order to avoid unnecessary memory allocation. Therefore, coordinates must be a tuple of 1d numpy arrays. This issue was already solved on point mass and prisms forward models, but the tesseroids forward model still stacks coordinates.

Are you willing to help implement and maintain this feature? Yes

Replace the Rio magnetic data with Great Britain open dataset

The Rio data is provided by CPRM but it doesn't really have a license. It would be best if we replaced it with an actually open dataset (with a permissive license).

The British Geological Survey provide a full aeromagnetic survey of Britain with a permissive open license. I would be happy to have a different dataset from outside North America and Europe but this might be the best we can get right now. We can always add another dataset later.

If this data file is too big (10s of Mb compressed), we can only host an interesting portion of it with nice anomalies.

What needs to be done:

  • Open and plot the dataset to see if it's worth using.
  • Convert the coordinates from OSGB36 to WGS84.
  • Keep only the useful columns (survey name/year, line number, longitude, latitude, flight height, IGRF90 corrected anomaly)
  • Save in compressed CVS and check the file size.
  • Commit archive and add to Pooch registry harmonica/datasets/registry.txt
  • Create a loading function and gallery example

Prisms forward modelling

Description of the desired feature

One of the most used geometries on three dimensional gravity modelling on projected coordinates are prisms.
Would be nice to port the prisms implementation from Fatiando 0.5 to Harmonica.

Are you willing to help implement and maintain this feature? Yes

Function to convert geocentric spherical coordinates back to geodetic

Description of the desired feature

In #28 we added the function geodetic_to_spherical to convert geodetic coordinates to geocentric spherical coordinates. It would be great to have the inverse transform as well (spherical_to_geodetic). The two functions should have similar docstrings/interfaces and include links to each other in their docstrings (using a See also section [example]).

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.