How to use harmonica - 10 common examples

region = [-5.5, -4.7, 57.8, 58.5]
inside = vd.inside((data.longitude, data.latitude), region)
data = data[inside]
print("Number of data points:", data.shape[0])
print("Mean height of observations:", data.altitude_m.mean())

# 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=1)

# 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 500 m spacing. The interpolation
# requires an extra coordinate (upward height). By passing in 1500 m, we're
# effectively upward-continuing the data (mean flight height is 500 m).
grid = eql.grid(spacing=500, data_names=["magnetic_anomaly"], extra_coords=1500)

# The grid is a xarray.Dataset with values, coordinates, and metadata
print("\nGenerated grid:\n", grid)

# Reduce region of the survey to speed things up
region = [-42.35, -42.10, -22.35, -22.15]
inside = vd.inside((data.longitude.values, data.latitude.values), region)
data = data[inside]
print("Number of data points:", data.shape[0])

# Project coordinates
projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())
data["easting"], data["northing"] = projection(
    data.longitude.values, data.latitude.values
)
coordinates = (data.easting, data.northing, data.altitude_m)

train, test = vd.train_test_split(coordinates, data.total_field_anomaly_nt, random_state=0)

eql = hm.EQLHarmonic(depth=1000, damping=10)
eql.fit(*train)
print("R² score on testing set:", eql.score(*test))

# Interpolate data into the regular grid at 200m above the sea level
grid = eql.grid(spacing=400, data_names=["magnetic_anomaly"], extra_coords=200)
# The grid is a xarray.Dataset with values, coordinates, and metadata
print(grid)

# Plot original magnetic anomaly
maxabs = vd.maxabs(data.total_field_anomaly_nt, grid.magnetic_anomaly.values)

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 6))
tmp = ax1.scatter(data.easting, data.northing, c=data.total_field_anomaly_nt, s=20,
                  vmin=-maxabs, vmax=maxabs, cmap="seismic")
plt.colorbar(tmp, ax=ax1, label="nT")
ax1.set_title("Observed Anomaly Magnetic data from Rio de Janeiro")

different heights. For a great number of applications we may need to interpolate these
data points into a regular grid. This can be done through an Equivalent Layer
interpolator. We will use :class:`harmonica.EQLHarmonic` to generate a set of point
sources beneath the observation points that predicts the observed data. Then these point
sources will be used to interpolate the data values into a regular grid at a constant
height.
"""
import matplotlib.pyplot as plt
import numpy as np
import pyproj
import verde as vd
import harmonica as hm


# Fetch magnetic anomaly data from Rio de Janeiro
data = hm.datasets.fetch_rio_magnetic()

# Reduce region of the survey to speed things up
region = [-42.35, -42.10, -22.35, -22.15]
inside = vd.inside((data.longitude.values, data.latitude.values), region)
data = data[inside]
print("Number of data points:", data.shape[0])

# Project coordinates
projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())
data["easting"], data["northing"] = projection(
    data.longitude.values, data.latitude.values
)
coordinates = (data.easting, data.northing, data.altitude_m)

train, test = vd.train_test_split(coordinates, data.total_field_anomaly_nt, random_state=0)

plot_include_source = True
plot_formats = ["png"]

# Sphinx project configuration
templates_path = ["_templates"]
exclude_patterns = ["_build", "**.ipynb_checkpoints"]
source_suffix = ".rst"
# The encoding of source files.
source_encoding = "utf-8-sig"
master_doc = "index"

# General information about the project
year = datetime.date.today().year
project = "Harmonica"
copyright = "2018-{}, The Harmonica Developers".format(year)
if len(full_version.split("+")) > 1 or full_version == "unknown":
    version = "dev"
else:
    version = full_version

# These enable substitutions using |variable| in the rst files
rst_epilog = """
.. |year| replace:: {year}
""".format(
    year=year
)

html_last_updated_fmt = "%b %d, %Y"
html_title = project
html_short_title = project
html_logo = "_static/harmonica-logo.png"
html_favicon = "_static/favicon.png"

# pylint: disable=missing-docstring,import-outside-toplevel
# Import functions/classes to make the public API
from . import version
from . import datasets
from . import synthetic
from .io import load_icgem_gdf
from .isostasy import isostasy_airy
from .gravity_corrections import bouguer_correction
from .forward.point_mass import point_mass_gravity
from .forward.tesseroid import tesseroid_gravity
from .forward.prism import prism_gravity
from .equivalent_layer.harmonic import EQLHarmonic, EQLHarmonicSpherical

# Get the version number through versioneer
__version__ = version.full_version


def test(doctest=True, verbose=True, coverage=False, figures=False):
    """
    Run the test suite.

    Uses `py.test `__ to discover and run the tests.

    Parameters
    ----------

    doctest : bool
        If ``True``, will run the doctests as well (code examples that start
        with a ``>>>`` in the docs).
    verbose : bool
        If ``True``, will print extra information during the test run.

blocked_mean = vd.BlockReduce(np.mean, spacing=0.2, drop_coords=False)
(longitude, latitude, elevation), gravity_data = blocked_mean.filter(
    (data.longitude, data.latitude, data.elevation), data.gravity,
)

# Compute gravity disturbance by removing the gravity of normal Earth
ellipsoid = bl.WGS84
gamma = ellipsoid.normal_gravity(latitude, height=elevation)
gravity_disturbance = gravity_data - gamma

# Convert data coordinates from geodetic (longitude, latitude, height) to
# spherical (longitude, spherical_latitude, radius).
coordinates = ellipsoid.geodetic_to_spherical(longitude, latitude, elevation)

# Create the equivalent layer
eql = hm.EQLHarmonicSpherical(damping=1e-3, relative_depth=10000)

# Fit the layer coefficients to the observed magnetic anomaly
eql.fit(coordinates, gravity_disturbance)

# 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, gravity_disturbance))

# Interpolate data on a regular grid with 0.2 degrees spacing. The
# interpolation requires an extra coordinate (radius). By passing in the
# maximum radius of the data, we're effectively upward-continuing the data.
# The grid will be defined in spherical coordinates.
grid = eql.grid(
    spacing=0.2, extra_coords=coordinates[-1].max(), data_names=["gravity_disturbance"],
)

===========

The geoid is the equipotential surface of the Earth's gravity potential that
coincides with mean sea level. It's often represented by "geoid heights", which
indicate the height of the geoid relative to the reference ellipsoid (WGS84 in
this case). Negative values indicate that the geoid is below the ellipsoid
surface and positive values that it is above. The data are on a regular grid
with 0.5 degree spacing and was generated from the spherical harmonic model
EIGEN-6C4 [Forste_etal2014]_.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import harmonica as hm

# Load the geoid grid
data = hm.datasets.fetch_geoid_earth()
print(data)

# Make a plot of data using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=100))
pc = data.geoid.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False)
plt.colorbar(
    pc, label="meters", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
)
ax.set_title("Geoid heights (EIGEN-6C4)")
ax.coastlines()
plt.tight_layout()
plt.show()

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 Great Britain
data = hm.datasets.fetch_britain_magnetic()

# Slice a smaller portion of the survey data to speed-up calculations for this
# example
region = [-5.5, -4.7, 57.8, 58.5]
inside = vd.inside((data.longitude, data.latitude), region)
data = data[inside]
print("Number of data points:", data.shape[0])
print("Mean height of observations:", data.altitude_m.mean())

# 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

(epsg:27700) coordinate system to WGS84 (epsg:4326) using to_crs function in
GeoPandas.

See the original data for more processing information.

If the file isn't already in your data directory, it will be downloaded
automatically.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import verde as vd
import harmonica as hm
import numpy as np

# Fetch the data in a pandas.DataFrame
data = hm.datasets.fetch_britain_magnetic()
print(data)

# Plot the observations in a Mercator map using Cartopy
fig = plt.figure(figsize=(7.5, 10))
ax = plt.axes(projection=ccrs.Mercator())
ax.set_title("Magnetic data from Great Britain", pad=25)
maxabs = np.percentile(data.total_field_anomaly_nt, 99)
tmp = ax.scatter(
    data.longitude,
    data.latitude,
    c=data.total_field_anomaly_nt,
    s=0.001,
    cmap="seismic",
    vmin=-maxabs,
    vmax=maxabs,
    transform=ccrs.PlateCarree(),

"""
Earth Gravity
=============

This is the magnitude of the gravity vector of the Earth (gravitational
+ centrifugal) at 10 km height. The data is on a regular grid with 0.5 degree
spacing at 10km ellipsoidal height. It was generated from the spherical
harmonic model EIGEN-6C4 [Forste_etal2014]_.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import harmonica as hm

# Load the gravity grid
data = hm.datasets.fetch_gravity_earth()
print(data)

# Make a plot of data using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=150))
pc = data.gravity.plot.pcolormesh(
    ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False
)
plt.colorbar(
    pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
)
ax.set_title("Gravity of the Earth (EIGEN-6C4)")
ax.coastlines()
plt.tight_layout()
plt.show()