diff --git a/examples/vector_uncoupled.py b/examples/vector_uncoupled.py index 343183026..acc6db163 100644 --- a/examples/vector_uncoupled.py +++ b/examples/vector_uncoupled.py @@ -1,13 +1,14 @@ """ -Gridding 2D vectors (uncoupled) -=============================== +Gridding 2D vectors +=================== We can use :class:`verde.Vector` to simultaneously process and grid all -components of vector data. Each component is processed and gridded separately (see -:class:`verde.VectorSpline2D` for a coupled alternative) but we have the convenience of -dealing with a single estimator. :class:`verde.Vector` can be combined with -:class:`verde.Trend`, :class:`verde.Spline`, and :class:`verde.Chain` to create a full -processing pipeline. +components of vector data. Each component is processed and gridded separately +(see `Erizo `__ for an elastically coupled +alternative) but we have the convenience of dealing with a single estimator. +:class:`verde.Vector` can be combined with :class:`verde.Trend`, +:class:`verde.Spline`, and :class:`verde.Chain` to create a full processing +pipeline. """ import matplotlib.pyplot as plt import cartopy.crs as ccrs diff --git a/examples/vectorspline2d.py b/examples/vectorspline2d.py deleted file mode 100644 index 6a84a4213..000000000 --- a/examples/vectorspline2d.py +++ /dev/null @@ -1,129 +0,0 @@ -""" -Gridding 2D vectors (coupled) -============================= - -One way of gridding vector data would be grid each component separately using -:class:`verde.Spline` and :class:`verde.Vector`. Alternatively, -:class:`verde.VectorSpline2D` can grid two components simultaneously in a way that -couples them through elastic deformation theory. This is particularly suited, though not -exclusive, to data that represent elastic/semi-elastic deformation, like horizontal GPS -velocities. -""" -import matplotlib.pyplot as plt -import cartopy.crs as ccrs -import numpy as np -import pyproj -import verde as vd - - -# Fetch the GPS data from the U.S. West coast. We'll grid only the horizontal components -# of the velocities -data = vd.datasets.fetch_california_gps() -coordinates = (data.longitude.values, data.latitude.values) -region = vd.get_region(coordinates) -# Use a Mercator projection because VectorSpline2D is a Cartesian gridder -projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean()) - -# Split the data into a training and testing set. We'll fit the gridder on the training -# set and use the testing set to evaluate how well the gridder is performing. -train, test = vd.train_test_split( - projection(*coordinates), (data.velocity_east, data.velocity_north), random_state=0 -) - -# We'll make a 15 arc-minute grid in the end. -spacing = 15 / 60 - -# Chain together a blocked mean to avoid aliasing, a polynomial trend to take care of -# the increase toward the coast, and finally the vector gridder using Poisson's ratio -# 0.5 to couple the two horizontal components. -chain = vd.Chain( - [ - ("mean", vd.BlockReduce(np.mean, spacing * 111e3)), - ("trend", vd.Vector([vd.Trend(degree=1) for i in range(2)])), - ("spline", vd.VectorSpline2D(poisson=0.5, mindist=10e3)), - ] -) -# Fit on the training data -chain.fit(*train) -# And score on the testing data. The best possible score is 1, meaning a perfect -# prediction of the test data. -score = chain.score(*test) -print("Cross-validation R^2 score: {:.2f}".format(score)) - -# Interpolate our horizontal GPS velocities onto a regular geographic grid and mask the -# data that are far from the observation points -grid_full = chain.grid( - region, spacing=spacing, projection=projection, dims=["latitude", "longitude"] -) -grid = vd.distance_mask( - (data.longitude, data.latitude), - maxdist=2 * spacing * 111e3, - grid=grid_full, - projection=projection, -) - -# Calculate residuals between the predictions and the original input data. Even though -# we aren't using regularization or regularly distributed forces, the prediction won't -# be perfect because of the BlockReduce operation. We fit the gridder on the reduced -# observations, not the original data. -predicted = chain.predict(projection(*coordinates)) -residuals = (data.velocity_east - predicted[0], data.velocity_north - predicted[1]) - -# Make maps of the original velocities, the gridded velocities, and the residuals -fig, axes = plt.subplots( - 1, 2, figsize=(12, 8), subplot_kw=dict(projection=ccrs.Mercator()) -) -crs = ccrs.PlateCarree() -# Plot the observed data and the residuals -ax = axes[0] -tmp = ax.quiver( - data.longitude.values, - data.latitude.values, - data.velocity_east.values, - data.velocity_north.values, - scale=0.3, - transform=crs, - width=0.001, - label="Velocities", -) -ax.quiverkey(tmp, 0.13, 0.18, 0.05, label="0.05 m/yr", coordinates="figure") -ax.quiver( - data.longitude.values, - data.latitude.values, - residuals[0].values, - residuals[1].values, - scale=0.3, - transform=crs, - color="r", - width=0.001, - label="Residuals", -) -ax.set_title("GPS horizontal velocities") -ax.legend(loc="lower left") -vd.datasets.setup_california_gps_map(ax) -# Plot the gridded data and the residuals -ax = axes[1] -tmp = ax.quiver( - grid.longitude.values, - grid.latitude.values, - grid.east_component.values, - grid.north_component.values, - scale=0.3, - transform=crs, - width=0.002, -) -ax.quiverkey(tmp, 0.63, 0.18, 0.05, label="0.05 m/yr", coordinates="figure") -ax.quiver( - data.longitude.values, - data.latitude.values, - residuals[0].values, - residuals[1].values, - scale=0.3, - transform=crs, - color="r", - width=0.001, -) -ax.set_title("Gridded velocities") -vd.datasets.setup_california_gps_map(ax) -plt.tight_layout() -plt.show() diff --git a/tutorials/trends.py b/tutorials/trends.py index 492d6f4eb..644c0bce3 100644 --- a/tutorials/trends.py +++ b/tutorials/trends.py @@ -2,11 +2,11 @@ Trend Estimation ================ -Trend estimation and removal is a common operation, particularly when dealing with -geophysical data. Moreover, some of the interpolation methods, like -:class:`verde.VectorSpline2D`, struggle with long-wavelength trends in the data. The -:class:`verde.Trend` class fits a 2D polynomial trend of arbitrary degree to the data -and can be used to remove it. +Trend estimation and removal is a common operation, particularly when dealing +with geophysical data. Moreover, some of the interpolation methods, like +:class:`verde.Spline`, can struggle with long-wavelength trends in the data. +The :class:`verde.Trend` class fits a 2D polynomial trend of arbitrary degree +to the data and can be used to remove it. """ import matplotlib.pyplot as plt import cartopy.crs as ccrs diff --git a/tutorials/vectors.py b/tutorials/vectors.py index d5ca8f2cd..63797fe61 100644 --- a/tutorials/vectors.py +++ b/tutorials/vectors.py @@ -281,59 +281,18 @@ plt.show() ######################################################################################## -# Another way of gridding 2-component vector data is using -# :class:`verde.VectorSpline2D`. This gridder uses linear elasticity theory to couple -# the two vector components. The degree of coupling can be controlled through the -# ``poisson`` parameter which sets the `Poisson's ratio -# `__ of the elastic medium. - -chain_coupled = vd.Chain( - [ - ("mean", vd.BlockMean(spacing=spacing * 111e3, uncertainty=True)), - ("trend", vd.Vector([vd.Trend(1), vd.Trend(1)])), - ("spline", vd.VectorSpline2D(poisson=0.5, damping=1e-10)), - ] -) -chain_coupled.fit(*train) -print(chain_coupled.score(*test)) - -######################################################################################## -# :class:`~verde.VectorSpline2D` generally gives better results when gridding GPS -# velocities, particularly for higher density grids and areas with sharp changes in -# velocity [SandwellWessel2016]_. Here, we won't see a big difference because of the -# low-density grid that we're making. - -grid_coupled = chain_coupled.grid( - region=region, - spacing=spacing, - projection=projection, - dims=["latitude", "longitude"], -) -grid_coupled = vd.distance_mask( - (data.longitude, data.latitude), - maxdist=spacing * 2 * 111e3, - grid=grid_coupled, - projection=projection, -) - -fig, axes = plt.subplots( - 1, 2, figsize=(9, 6.5), subplot_kw=dict(projection=ccrs.Mercator()) -) -crs = ccrs.PlateCarree() -titles = ["Gridded velocity (uncoupled)", "Gridded velocity (coupled)"] -grids = [grid, grid_coupled] -for ax, grd, title in zip(axes, grids, titles): - ax.set_title(title) - tmp = ax.quiver( - grd.longitude.values, - grd.latitude.values, - grd.east_component.values, - grd.north_component.values, - scale=0.3, - transform=crs, - width=0.002, - ) - vd.datasets.setup_california_gps_map(ax) -ax.quiverkey(tmp, 0.15, 0.15, 0.05, label="0.05 m/yr", coordinates="figure") -plt.tight_layout() -plt.show() +# GPS/GNSS data +# +++++++++++++ +# +# For some types of vector data, like GPS/GNSS displacements, the vector +# components are coupled through elasticity. In these cases, elastic Green's +# functions can be used to achieve better interpolation results. The `Erizo +# package `__ implements some of these +# Green's functions. +# +# .. warning:: +# +# The :class:`verde.VectorSpline2D` class implemented an elastically +# coupled Green's function but it is deprecated and will be removed in +# Verde v2.0.0. Please use the implementation in the `Erizo +# `__ package instead. diff --git a/verde/vector.py b/verde/vector.py index db4f8049e..ff08fad36 100644 --- a/verde/vector.py +++ b/verde/vector.py @@ -1,6 +1,8 @@ """ Classes for dealing with vector data. """ +import warnings + import numpy as np from sklearn.utils.validation import check_is_fitted @@ -17,6 +19,9 @@ from .utils import dummy_jit as jit +# Otherwise, DeprecationWarning won't be shown, kind of defeating the purpose. +warnings.simplefilter("default") + # Default arguments for numba.jit JIT_ARGS = dict(nopython=True, target="cpu", fastmath=True, parallel=True) @@ -141,6 +146,13 @@ class VectorSpline2D(BaseGridder): r""" Elastically coupled interpolation of 2-component vector data. + .. warning:: + + The :class:`~verde.VectorSpline2D` class is deprecated and will be + removed in Verde v2.0.0. Its usage is restricted to GPS/GNSS data and + not in the general scope of Verde. Please use the implementation in the + `Erizo `__ package instead. + This gridder assumes Cartesian coordinates. Uses the Green's functions based on elastic deformation from @@ -221,6 +233,12 @@ def __init__( self.damping = damping self.force_coords = force_coords self.engine = engine + warnings.warn( + "VectorSpline2D is deprecated and will be removed in Verde v2.0.0." + " Please use the implementation in the Erizo package instead " + "(https://github.com/fatiando/erizo).", + DeprecationWarning, + ) def fit(self, coordinates, data, weights=None): """