Code Monkey home page Code Monkey logo

infresnel's Introduction

infresnel

Documentation Status Build Status Code style: black Imports: isort Binder
DOI

infresnel is a Python package which facilitates computation of the Fresnel number for infrasound applications — specifically, propagation over topography at local to regional scales.

Background and motivation

The dimensionless Fresnel number $N$ for an acoustic wave with wavelength $\lambda$ is given by $$N = \frac{(R_\mathrm{f} - R_\mathrm{d})}{\lambda / 2},$$ where $R_\mathrm{f}$ is the length of the shortest diffracted path over topography, and $R_\mathrm{d}$ the length of the direct path (i.e., line-of-sight slant distance), from source to receiver (Maher et al., 2021). In words, the Fresnel number is the extra distance a wave must travel due to topography, normalized by half a wavelength. The path length difference $(R_\mathrm{f} - R_\mathrm{d})$ can additionally be used to estimate travel time delays due to diffraction over topography. The travel time delay $\delta_\mathrm{f}$ is given by $$\delta_\mathrm{f} = \frac{(R_\mathrm{f} - R_\mathrm{d})}{c},$$ where $c$ is the estimated acoustic wavespeed assuming a homogenous atmosphere. $\delta_\mathrm{f}$ has been shown to be comparable to travel time differences computed using finite-difference time-domain simulations — e.g., see Fig. 2B in Fee et al. (2021).

$N$ can be used to estimate the loss of infrasound energy due to diffraction over topography — the "insertion loss" — via application of empirical scaling relationships. See Maher et al. (2021) for a thorough discussion of these various empirical relationships, and their limitations, in the context of volcano infrasound.

These are simple equations, but the practical computation of the quantity $(R_\mathrm{f} - R_\mathrm{d})$ is somewhat involved. The goal of infresnel is to make this computation as quick and convenient as possible.

Quickstart

  1. Obtain

    git clone https://github.com/liamtoney/infresnel.git
    cd infresnel
    
  2. Create environment, install, and activate (install conda first, if necessary)

    conda env create
    conda activate infresnel
    
  3. Run using the Python interpreter

    python
    >>> from infresnel import calculate_paths

Usage

Interactive API reference documentation for infresnel can be found here.

Additionally, several interactive notebooks containing usage examples are included in notebooks. To open the notebooks, with your conda environment activated run

jupyter lab notebooks

Alternatively, you can run these notebooks in your browser — without installing infresnel — by navigating to the project's Binder.

Assumptions

infresnel calculates path length differences using elevation profiles. This means that all diffraction is assumed to take place in the vertical plane between source and receiver. One can easily construct scenarios where this assumption is violated. For example, consider a column of rock much taller than it is wide, located directly between source and receiver. For this scenario, infresnel would predict a large path length difference for waves traveling over the top of the column — but in reality, wavefronts diffract laterally around the column. The true travel time from source to receiver is thus much smaller than what infresnel predicts for this scenario.

Citing

If you use infresnel in research that leads to a published manuscript, please cite the tool:

Toney, L. (2024). infresnel (v0.3.0). Zenodo. https://doi.org/10.5281/zenodo.11176356

To cite a previous version of infresnel, go to the Zenodo page, select the version, and use the "Cite as" tool there.

Installation details

infresnel's dependencies are:

  • colorcet (for perceptually accurate colormaps to use in the creation of GeoTIFF path length difference grids)
  • ipympl
  • ipywidgets
  • joblib (for parallel processing)
  • JupyterLab (for running the interactive .ipynb notebooks)
  • Matplotlib (for applying colormaps to GeoTIFFs and for generally plotting results)
  • Numba (for computational acceleration)
  • pip
  • PyGMT (for simplified SRTM data downloading and caching)
  • rioxarray (for DEM file I/O, reprojection, and elevation profile interpolation)
  • tqdm (for displaying progress bars)

These dependencies are listed in environment.yml, and they are installed in step 2 above.

You might want to install infresnel into an existing conda environment, instead of making a new one. In this case, after step 1 above run

conda activate <existing_environment>
pip install --editable .

which uses pip to install infresnel's dependencies, if you don't already have them installed in your existing environment.

In either case, your installation will be "editable." This means that you can modify the source code in your local infresnel directory — or run a git pull to update with any new remote changes — and the installed package will be updated. To instead use a specific release of infresnel, run

git switch -d vX.Y.Z

between steps 1 and 2, where vX.Y.Z is the release version (e.g., v0.1.0). You can switch back with git switch main.

Contributing

If you notice a bug with infresnel (or if you'd like to request/propose a new feature), please create an issue on GitHub (preferred) or email me at [email protected]. You are also welcome to create a pull request.

To install infresnel's development packages, with your conda environment activated run

pip install --requirement requirements.txt
nbdime config-git --enable  # Configure Jupyter Notebook diffs

References

Fee, D., Toney, L., Kim, K., Sanderson, R. W., Iezzi, A. M., Matoza, R. S., De Angelis, S., Jolly, A. D., Lyons, J. J., & Haney, M. (2021). Local explosion detection and infrasound localization by reverse time migration using 3-D finite-difference wave propagation. Frontiers in Earth Science, 9. https://doi.org/10.3389/feart.2021.620813

Maher, S. P., Matoza, R. S., de Groot-Hedlin, C., Kim, K., & Gee, K. L. (2021). Evaluating the applicability of a screen diffraction approximation to local volcano infrasound. Volcanica, 4(1), 67–85. https://doi.org/10.30909/vol.04.01.6785

Acknowledgements

This work was supported by the Nuclear Arms Control Technology (NACT) program at the Defense Threat Reduction Agency (DTRA). Cleared for release.

infresnel's People

Contributors

liamtoney avatar

Stargazers

 avatar  avatar  avatar  avatar

Watchers

 avatar

infresnel's Issues

Incorrect UTM CRS

I get an incorrect UTM geodetic(?) CRS for a scenario in Italy. It seems to come from

utm_crs = _estimate_utm_crs(src_lat, src_lon, datum_name='WGS 84')

where

utm_crs = _estimate_utm_crs(15.2122, 38.7946, datum_name='WGS 84')

which should return a CRS of 33 but instead returns one of 37.

DEMs with NaN values produce "cliffs" that artificially blow up path length difference values

Basically, we should implement a "safety buffer" around the valid-elevation-data borders of such DEMs. I.e., we don't want to compute paths for points near that border. This is mainly relevant for calculate_paths_grid().

Not a hugely critical issue since it's usually pretty obvious when this happens. Applying e.g. vmax=np.nanpercentile(path_length_differences, 99) when creating an image plot works pretty well to suppress these artificially high outliers.

Perhaps just editing the tolerance in this part of the code is enough? E.g. increasing to 2 * mean_resolution or something... easy to try this out

ZeroDivisionError when grid spacing < profile spacing

Sometimes (usually when running infresnel to produce grids) a ZeroDivisionError occurs. The error is raised during this line:

return (z[-1] - z[0]) / (d[-1] - d[0]) * (d - d[0]) + z[0] # d - d[0] for intercept

Evidently, sometimes the quantity d[-1] - d[0] is 0. This can occur when the profile spacing (which is currently automatically determined from the DEM spacing) is larger than the requested grid spacing.

Optimize code

The main computational expense involved in the code is the DEM interpolation step. There may be ways to expedite this that don't take too much effort to implement.

Handling NaN values in input DEMs

Not a Number (NaN; e.g., np.nan) values can make their way into input DEMs for infresnel in a few ways. There are no SRTM data at high latitudes (beyond 60° N and 56° S), so the automatic DEM downloader may return NaN values (check this!). Also, a user-supplied DEM may have NaN values either because it contains holes (yikes) or because it has a non-rectangular geometry in the UTM projection (more likely; see the Yasur DEM for example). I see three cases to address:

  1. The DEM is all NaN values. This might occur if dem_file=None and we are totally beyond the SRTM latitude range. We should error out here.
  2. The DEM has no NaN values. This is generally the case for SRTM data totally within the latitude range, since we are clipping from a much larger DEM. It can also be the case for well-trimmed user-supplied DEMs.
  3. The DEM is partially composed of NaN values. This can be the case for SRTM data on the "boundary" of the latitude range, or for user-supplied DEMs with holes or non-rectangular geometries as mentioned above. We should warn the user here, since NaN values in the DEM pose problems. We can print the percentage of the DEM that is NaN as a useful warning...

For case (3) above, even without holes in the DEM we still have the challenge of fitting the RectBivariateSpline, which cannot have any NaN values in the input z. Currently we just replace NaN with 0 everywhere to get around this. However, this is not the best idea since it produces large discontinuities ("cliffs") in elevation where the spline fitting becomes extremely unstable. This can cause the evaluated elevation profiles to have artifacts — which bleeds into the path length computations.

Even though we can skip computing paths for receivers located in NaN regions, receivers just on the edge of the valid (non-NaN) portion of the DEM can still experience edge effects. There are two approaches we can take to avoid this behavior, both which only apply for case (3):

  • Instead of filling with 0, we can fill with a more intelligent value such as the median of the DEM data. This will still result in discontinuities but may be sufficient to suppress edge effects in most cases. This is fast.
  • We can perform a nearest-neighbor extrapolation, via e.g. scipy.interpolate.griddata(). This basically matches the area around the valid portion with similar values. This is obviously not valid far from the valid region, but it should be very effective in stabilizing the spline fitting process. This is slow. It might make sense to downsample the evaluation points xi for speed, then upsample them to the DEM resolution before adding as the "background" to the DEM to fill the NaN values. This is significantly more work and adds complexity to the code.

Tentative action items, based on the above...

  • Add a check for NaN values in the DEM, to determine which of the 3 cases we have.
  • If case (3), print the percentage of NaN values as a warning, then apply the first of the two approaches outlined above prior to fitting the spline.

We can see if this is successful to suppress artifacts; if not then the second of the two approaches above may need to be implemented.

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.