forked from fatiando/harmonica
-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #1 from fatiando/master
bring recent changes from fatiando/harmonica master
- Loading branch information
Showing
28 changed files
with
1,469 additions
and
156 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,41 @@ | ||
| """ | ||
| Land Gravity Data from South Africa | ||
| =================================== | ||
| Land gravity survey performed in January 1986 within the boundaries of the Republic of | ||
| South Africa. The data was made available by the `National Centers for Environmental | ||
| Information (NCEI) <https://www.ngdc.noaa.gov/>`__ (formerly NGDC) and are in the | ||
| `public domain <https://www.ngdc.noaa.gov/ngdcinfo/privacy.html#copyright-notice>`__. | ||
| The entire dataset is stored in a :class:`pandas.DataFrame` with columns: longitude, | ||
| latitude, elevation (above sea level) and gravity(mGal). See the documentation for | ||
| :func:`harmonica.datasets.fetch_south_africa_gravity` for more information. | ||
| """ | ||
| import matplotlib.pyplot as plt | ||
| import cartopy.crs as ccrs | ||
| import verde as vd | ||
| import harmonica as hm | ||
|
|
||
| # Fetch the data in a pandas.DataFrame | ||
| data = hm.datasets.fetch_south_africa_gravity() | ||
| print(data) | ||
|
|
||
| # Plot the observations in a Mercator map using Cartopy | ||
| fig = plt.figure(figsize=(6.5, 5)) | ||
| ax = plt.axes(projection=ccrs.Mercator()) | ||
| ax.set_title("Observed gravity data from South Africa", pad=25) | ||
| tmp = ax.scatter( | ||
| data.longitude, | ||
| data.latitude, | ||
| c=data.gravity, | ||
| s=0.8, | ||
| cmap="viridis", | ||
| transform=ccrs.PlateCarree(), | ||
| ) | ||
| plt.colorbar( | ||
| tmp, ax=ax, label="observed gravity [mGal]", aspect=50, pad=0.1, shrink=0.92 | ||
| ) | ||
| ax.set_extent(vd.get_region((data.longitude, data.latitude))) | ||
| ax.gridlines(draw_labels=True) | ||
| ax.coastlines() | ||
| plt.tight_layout() | ||
| plt.show() |
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,2 @@ | ||
| Equivalent Layers | ||
| ----------------- |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,97 @@ | ||
| """ | ||
| Gridding and upward continuation | ||
| ================================ | ||
| Most potential field surveys gather data along scattered and uneven flight lines or | ||
| ground measurements. For a great number of applications we may need to interpolate these | ||
| data points onto a regular grid at a constant altitude. Upward-continuation is also a | ||
| routine task for smoothing, noise attenuation, source separation, etc. | ||
| Both tasks can be done simultaneously through an *equivalent layer* [Dampney1969]_. We | ||
| will use :class:`harmonica.EQLHarmonic` to estimate the coefficients of a set of point | ||
| sources (the equivalent layer) that fit the observed data. The fitted layer can then be | ||
| used to predict data values wherever we want, like on a grid at a certain altitude. The | ||
| sources for :class:`~harmonica.EQLHarmonic` in particular are placed one beneath each | ||
| data point at a constant relative depth from the elevation of the data point following | ||
| [Cooper2000]_. | ||
| The advantage of using an equivalent layer is that it takes into account the 3D nature | ||
| of the observations, not just their horizontal positions. It also allows data | ||
| uncertainty to be taken into account and noise to be suppressed though the least-squares | ||
| fitting process. The main disadvantage is the increased computational load (both in | ||
| terms of time and memory). | ||
| """ | ||
| import matplotlib.pyplot as plt | ||
| import numpy as np | ||
| import pyproj | ||
| import verde as vd | ||
| import harmonica as hm | ||
|
|
||
|
|
||
| # Fetch the sample total-field magnetic anomaly data from Rio de Janeiro | ||
| data = hm.datasets.fetch_rio_magnetic() | ||
|
|
||
| # Slice a smaller portion of the survey data to speed-up calculations for this example | ||
| region = [-42.35, -42.10, -22.35, -22.15] | ||
| inside = vd.inside((data.longitude, data.latitude), region) | ||
| data = data[inside] | ||
| print("Number of data points:", data.shape[0]) | ||
|
|
||
| # Since this is a small area, we'll project our data and use Cartesian coordinates | ||
| projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean()) | ||
| easting, northing = projection(data.longitude.values, data.latitude.values) | ||
| coordinates = (easting, northing, data.altitude_m) | ||
|
|
||
| # Create the equivalent layer. We'll use the default point source configuration at a | ||
| # constant relative depth beneath each observation point. The damping parameter helps | ||
| # smooth the predicted data and ensure stability. | ||
| eql = hm.EQLHarmonic(relative_depth=1000, damping=10) | ||
|
|
||
| # Fit the layer coefficients to the observed magnetic anomaly. | ||
| eql.fit(coordinates, data.total_field_anomaly_nt) | ||
|
|
||
| # Evaluate the data fit by calculating an R² score against the observed data. | ||
| # This is a measure of how well layer the fits the data NOT how good the interpolation | ||
| # will be. | ||
| print("R² score:", eql.score(coordinates, data.total_field_anomaly_nt)) | ||
|
|
||
| # Interpolate data on a regular grid with 400 m spacing. The interpolation requires an | ||
| # extra coordinate (upward height). By passing in 500 m, we're effectively | ||
| # upward-continuing the data (mean flight height is 200 m). | ||
| grid = eql.grid(spacing=400, data_names=["magnetic_anomaly"], extra_coords=500) | ||
|
|
||
| # The grid is a xarray.Dataset with values, coordinates, and metadata | ||
| print("\nGenerated grid:\n", grid) | ||
|
|
||
| # Plot original magnetic anomaly and the gridded and upward-continued version | ||
| fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 6), sharey=True) | ||
|
|
||
| # Get the maximum absolute value between the original and gridded data so we can use the | ||
| # same color scale for both plots and have 0 centered at the white color. | ||
| maxabs = vd.maxabs(data.total_field_anomaly_nt, grid.magnetic_anomaly.values) | ||
|
|
||
| ax1.set_title("Observed magnetic anomaly data") | ||
| tmp = ax1.scatter( | ||
| easting, | ||
| northing, | ||
| c=data.total_field_anomaly_nt, | ||
| s=20, | ||
| vmin=-maxabs, | ||
| vmax=maxabs, | ||
| cmap="seismic", | ||
| ) | ||
| plt.colorbar(tmp, ax=ax1, label="nT", pad=0.05, aspect=40, orientation="horizontal") | ||
|
|
||
| ax2.set_title("Gridded and upward-continued") | ||
| tmp = grid.magnetic_anomaly.plot.pcolormesh( | ||
| ax=ax2, | ||
| add_colorbar=False, | ||
| add_labels=False, | ||
| vmin=-maxabs, | ||
| vmax=maxabs, | ||
| cmap="seismic", | ||
| ) | ||
| plt.colorbar(tmp, ax=ax2, label="nT", pad=0.05, aspect=40, orientation="horizontal") | ||
|
|
||
| plt.tight_layout(pad=0) | ||
| plt.show() |
Oops, something went wrong.