Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion examples/equivalent_sources/block_averaged_sources.py
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,8 @@
# requires the height of the grid points (upward coordinate). By passing in
# 1500 m, we're effectively upward-continuing the data (mean flight height is
# 500 m).
grid = eqs.grid(upward=1500, spacing=500, data_names=["magnetic_anomaly"])
grid_coords = vd.grid_coordinates(region=region, spacing=500, extra_coords=1500)
grid = eqs.grid(coordinates=grid_coords, data_names=["magnetic_anomaly"])

# The grid is a xarray.Dataset with values, coordinates, and metadata
print("\nGenerated grid:\n", grid)
Expand Down
3 changes: 2 additions & 1 deletion examples/equivalent_sources/cartesian.py
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,8 @@
# requires the height of the grid points (upward coordinate). By passing in
# 1500 m, we're effectively upward-continuing the data (mean flight height is
# 500 m).
grid = eqs.grid(upward=1500, spacing=500, data_names=["magnetic_anomaly"])
grid_coords = vd.grid_coordinates(region=region, spacing=500, extra_coords=1500)
grid = eqs.grid(coordinates=grid_coords, data_names=["magnetic_anomaly"])

# The grid is a xarray.Dataset with values, coordinates, and metadata
print("\nGenerated grid:\n", grid)
Expand Down
7 changes: 2 additions & 5 deletions examples/equivalent_sources/gradient_boosted.py
Original file line number Diff line number Diff line change
Expand Up @@ -91,11 +91,8 @@
# Interpolate data on a regular grid with 2 km spacing. The interpolation
# requires the height of the grid points (upward coordinate). By passing in
# 1000 m, we're effectively upward-continuing the data.
grid = eqs_gb.grid(
upward=1000,
spacing=2e3,
data_names="gravity_disturbance",
)
grid_coords = vd.grid_coordinates(region=region, spacing=2e3, extra_coords=1000)
grid = eqs_gb.grid(coordinates=grid_coords, data_names="gravity_disturbance")
print(grid)

# Plot the original gravity disturbance and the gridded and upward-continued
Expand Down
12 changes: 4 additions & 8 deletions examples/equivalent_sources/spherical.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,11 +62,10 @@
# passing in the maximum radius of the data, we're effectively
# upward-continuing the data. The grid will be defined in spherical
# coordinates.
grid = eqs.grid(
upward=coordinates[-1].max(),
spacing=0.2,
data_names=["gravity_disturbance"],
)
region = vd.get_region(coordinates) # get the region boundaries
upward = coordinates[-1].max()
grid_coords = vd.grid_coordinates(region=region, spacing=0.2, extra_coords=upward)
grid = eqs.grid(coordinates=grid_coords, data_names=["gravity_disturbance"])

# The grid is a xarray.Dataset with values, coordinates, and metadata
print("\nGenerated grid:\n", grid)
Expand All @@ -79,9 +78,6 @@
# color.
maxabs = vd.maxabs(gravity_disturbance, grid.gravity_disturbance.values)

# Get the region boundaries
region = vd.get_region(coordinates)

# Plot observed and gridded gravity disturbance
fig, (ax1, ax2) = plt.subplots(
nrows=1,
Expand Down
78 changes: 47 additions & 31 deletions harmonica/equivalent_sources/cartesian.py
Original file line number Diff line number Diff line change
Expand Up @@ -377,10 +377,7 @@ def jacobian(self, coordinates, points):

def grid(
self,
upward,
region=None,
shape=None,
spacing=None,
coordinates,
dims=None,
data_names=None,
projection=None,
Expand All @@ -389,31 +386,27 @@ def grid(
"""
Interpolate the data onto a regular grid.

The grid can be specified by either the number of points in each
dimension (the *shape*) or by the grid node spacing. See
:func:`verde.grid_coordinates` for details. All grid points will be
located at the same `upward` coordinate. Other arguments for
:func:`verde.grid_coordinates` can be passed as extra keyword arguments
(``kwargs``) to this method.
The coordinates of the regular grid must be passed through the
``coordinates`` argument as a tuple containing three arrays in the
following order: ``(easting, nothing, upward)``. They can be easily
created through the :func:`verde.grid_coordinates` function. If the
grid points must be all at the same height, it can be specified in the
``extra_coords`` argument of :func:`verde.grid_coordinates`.

If the interpolator collected the input data region, then it will be
used if ``region=None``. Otherwise, you must specify the grid region.
Use the *dims* and *data_names* arguments to set custom names for the
dimensions and the data field(s) in the output :class:`xarray.Dataset`.
Default names will be provided if none are given.

Parameters
----------
upward : float
Upward coordinate of the grid points.
region : list = [W, E, S, N]
The west, east, south, and north boundaries of a given region.
shape : tuple = (n_north, n_east) or None
The number of points in the South-North and West-East directions,
respectively.
spacing : tuple = (s_north, s_east) or None
The grid spacing in the South-North and West-East directions,
respectively.
coordinates : tuple of arrays
Tuple of arrays containing the coordinates of the grid in the
following order: (easting, northing, upward).
The easting and northing arrays could be 1d or 2d arrays, if
they are 2d they must be part of a meshgrid.
The upward array should be a 2d array with the same shape of
easting and northing (if they are 2d arrays) or with a shape of
``(northing.size, easting.size)`` (if they are 1d arrays).
dims : list or None
The names of the northing and easting data dimensions,
respectively, in the output grid. Default is determined from the
Expand All @@ -439,22 +432,45 @@ def grid(
The interpolated grid. Metadata about the interpolator is written
to the ``attrs`` attribute.

See also
--------
:func:`verde.grid_coordinates`

"""
# We override the grid method from BaseGridder so it takes the upward
# coordinate as a positional argument.
# We override the grid method from BaseGridder to change the docstring
# and to make it work only with the `coordinates` argument (no region,
# shape or spacing)

# Raise ValueError if any deprecated argument has been passed
deprecated_args = (
"upward" in kwargs,
"shape" in kwargs,
"region" in kwargs,
"spacing" in kwargs,
)
if any(deprecated_args):
raise ValueError(
"The 'upward', 'region', 'shape' and 'spacing' arguments have been "
+ "deprecated. "
+ "Please pass the coordinates of the target grid through the "
+ "'coordinates' argument."
)

# Raise warning if any kwargs has been passed
if kwargs:
args = "'" + "', '".join(list(kwargs.keys())) + "'"
warnings.warn(
f"The {args} arguments are being ignored. The 'grid' method "
+ "will not take any keyword arguments in the next Harmonica release",
FutureWarning,
)

# Ignore extra_coords if passed
pop_extra_coords(kwargs)
# Grid data
grid = super().grid(
region=region,
shape=shape,
spacing=spacing,
coordinates=coordinates,
dims=dims,
data_names=data_names,
projection=projection,
extra_coords=upward,
**kwargs,
)
return grid

Expand Down
83 changes: 47 additions & 36 deletions harmonica/equivalent_sources/spherical.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@
from .utils import (
jacobian_numba_parallel,
jacobian_numba_serial,
pop_extra_coords,
predict_numba_parallel,
predict_numba_serial,
)
Expand Down Expand Up @@ -239,49 +238,42 @@ def jacobian(self, coordinates, points, dtype="float64"):

def grid(
self,
upward,
region=None,
shape=None,
spacing=None,
coordinates,
dims=None,
data_names=None,
**kwargs,
):
"""
Interpolate the data onto a regular grid.

The grid can be specified by either the number of points in each
dimension (the *shape*) or by the grid node spacing. See
:func:`verde.grid_coordinates` for details. All grid points will be
located at the same `upward` coordinate. Other arguments for
:func:`verde.grid_coordinates` can be passed as extra keyword arguments
(``kwargs``) to this method.
The coordinates of the regular grid must be passed through the
``coordinates`` argument as a tuple containing three arrays in the
following order: ``(longitude, latitude, radius)``. They can be easily
created through the :func:`verde.grid_coordinates` function. If the
grid points must be all at the same radius, it can be specified in the
``extra_coords`` argument of :func:`verde.grid_coordinates`.

If the interpolator collected the input data region, then it will be
used if ``region=None``. Otherwise, you must specify the grid region.
Use the *dims* and *data_names* arguments to set custom names for the
dimensions and the data field(s) in the output :class:`xarray.Dataset`.
Default names will be provided if none are given.

Parameters
----------
upward : float
Upward coordinate of the grid points.
region : list = [W, E, S, N]
The west, east, south, and north boundaries of a given region.
shape : tuple = (n_north, n_east) or None
The number of points in the South-North and West-East directions,
respectively.
spacing : tuple = (s_north, s_east) or None
The grid spacing in the South-North and West-East directions,
respectively.
coordinates : tuple of arrays
Tuple of arrays containing the coordinates of the grid in the
following order: (longitude, latitude, radius).
The longitude and latitude arrays could be 1d or 2d arrays, if
they are 2d they must be part of a meshgrid.
The radius array should be a 2d array with the same shape of
longitude and latitude (if they are 2d arrays) or with a shape of
``(latitude.size, longitude.size)`` (if they are 1d arrays).
dims : list or None
The names of the northing and easting data dimensions,
The names of the latitude and longitude data dimensions,
respectively, in the output grid. Default is determined from the
``dims`` attribute of the class. Must be defined in the following
order: northing dimension, easting dimension.
**NOTE: This is an exception to the "easting" then
"northing" pattern but is required for compatibility with xarray.**
order: latitude dimension, longitude dimension.
**NOTE: This is an exception to the "longitude" then
"latitude" pattern but is required for compatibility with xarray.**
data_names : list of None
The name(s) of the data variables in the output grid. Defaults to
``['scalars']``.
Expand All @@ -293,23 +285,42 @@ def grid(
to the ``attrs`` attribute.

"""
# We override the grid method from BaseGridder so it takes the upward
# coordinate as a positional argument.
# We override the grid method from BaseGridder to change the docstring
# and to make it work only with the `coordinates` argument (no region,
# shape or spacing)

# Raise ValueError if any deprecated argument has been passed
deprecated_args = (
"upward" in kwargs,
"shape" in kwargs,
"region" in kwargs,
"spacing" in kwargs,
)
if any(deprecated_args):
raise ValueError(
"The 'upward', 'region', 'shape' and 'spacing' arguments have been "
+ "deprecated. "
+ "Please pass the coordinates of the target grid through the "
+ "'coordinates' argument."
)

# Raise warning if any kwargs has been passed
if kwargs:
args = "'" + "', '".join(list(kwargs.keys())) + "'"
warnings.warn(
f"The {args} arguments are being ignored. The 'grid' method "
+ "will not take any keyword arguments in the next Harmonica release",
FutureWarning,
)

# Ignore extra_coords if passed
pop_extra_coords(kwargs)
# Grid data
# We always pass projection=None because that argument it's intended to
# be used only with Cartesian gridders.
grid = super().grid(
region=region,
shape=shape,
spacing=spacing,
coordinates=coordinates,
dims=dims,
data_names=data_names,
projection=None,
extra_coords=upward,
**kwargs,
)
return grid

Expand Down
Loading