How to use the verde.datasets.fetch_california_gps function in verde

: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

Verde provides the :class:`verde.Trend` class to estimate a polynomial trend and the
:class:`verde.Vector` class to apply any combination of estimators to each component
of vector data, like GPS velocities. You can access each component as a separate
(fitted) :class:`verde.Trend` instance or operate on all vector components directly
using using :meth:`verde.Vector.predict`, :meth:`verde.Vector.grid`, etc, or
chaining it with a vector interpolator using :class:`verde.Chain`.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import verde as vd


# Fetch the GPS data from the U.S. West coast. The data has a strong trend toward the
# North-West because of the relative movement along the San Andreas Fault System.
data = vd.datasets.fetch_california_gps()

# We'll fit a degree 2 trend on both the East and North components and weight the data
# using the inverse of the variance of each component.
# Note: Never use [Trend(...)]*2 as an argument to Vector. This creates references
# to the same Trend instance and will mess up the fitting.
trend = vd.Vector([vd.Trend(degree=2) for i in range(2)])
weights = vd.variance_to_weights((data.std_east ** 2, data.std_north ** 2))
trend.fit(
    coordinates=(data.longitude, data.latitude),
    data=(data.velocity_east, data.velocity_north),
    weights=weights,
)
print("Vector trend estimator:", trend)

# The separate Trend objects for each component can be accessed through the 'components'
# attribute. You could grid them individually if you wanted.

Verde provides the :class:`verde.Trend` class to estimate a polynomial trend and the
:class:`verde.Vector` class to apply any combination of estimators to each component
of vector data, like GPS velocities. You can access each component as a separate
(fitted) :class:`verde.Trend` instance or operate on all vector components directly
using using :meth:`verde.Vector.predict`, :meth:`verde.Vector.grid`, etc, or
chaining it with a vector interpolator using :class:`verde.Chain`.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import verde as vd


# Fetch the GPS data from the U.S. West coast. The data has a strong trend toward the
# North-West because of the relative movement along the San Andreas Fault System.
data = vd.datasets.fetch_california_gps()

# We'll fit a degree 2 trend on both the East and North components and weight the data
# using the inverse of the variance of each component.
# Note: Never use [Trend(...)]*2 as an argument to Vector. This creates references
# to the same Trend instance and will mess up the fitting.
trend = vd.Vector([vd.Trend(degree=2) for i in range(2)])
weights = vd.variance_to_weights((data.std_east ** 2, data.std_north ** 2))
trend.fit(
    coordinates=(data.longitude, data.latitude),
    data=(data.velocity_east, data.velocity_north),
    weights=weights,
)
print("Vector trend estimator:", trend)

# The separate Trend objects for each component can be accessed through the 'components'
# attribute. You could grid them individually if you wanted.

==================================

Sometimes data has outliers or less reliable points that might skew a blocked mean or
even a median. If the reduction function can take a ``weights`` argument, like
``numpy.average``, you can pass in weights to :class:`verde.BlockReduce` to lower the
influence of the offending data points. However, :class:`verde.BlockReduce` can't
produce weights for the blocked data (for use by a gridder, for example). If you want to
produced blocked weights as well, use :class:`verde.BlockMean`.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import verde as vd

# We'll test this on the California vertical GPS velocity data
data = vd.datasets.fetch_california_gps()

# We'll add some random extreme outliers to the data
outliers = np.random.RandomState(2).randint(0, data.shape[0], size=20)
data.velocity_up[outliers] += 0.08
print("Index of outliers:", outliers)

# Create an array of weights and set the weights for the outliers to a very low value
weights = np.ones_like(data.velocity_up)
weights[outliers] = 1e-5

# Now we can block average the points with and without weights to compare the outputs.
reducer = vd.BlockReduce(reduction=np.average, spacing=30 / 60, center_coordinates=True)
coordinates, no_weights = reducer.filter(
    (data.longitude, data.latitude), data.velocity_up
)
__, with_weights = reducer.filter(

low weight for data that have large uncertainties but don't vary much inside the block.

The weighted variance should be used when the data vary a lot in each block (high
variance) but have very similar uncertainties. The variance will be large when there is
a lot of variability in the data that isn't due to the uncertainties. This is also the
best choice if your data weights aren't ``1/uncertainty**2``.
"""
import matplotlib.pyplot as plt
from matplotlib.colors import PowerNorm, LogNorm
import cartopy.crs as ccrs
import numpy as np
import verde as vd

# We'll test this on the California vertical GPS velocity data because it comes with the
# uncertainties
data = vd.datasets.fetch_california_gps()
coordinates = (data.longitude, data.latitude)

# We'll calculate the mean on large blocks to show the effect of the different weighting
# schemes
spacing = 30 / 60
# It's important that the weights are given as 1/sigma**2 for the uncertainty
# propagation. In this case, you should not use verde.variance_to_weights because it
# would normalize the weights.
weights = 1 / data.std_up ** 2
reducer = vd.BlockMean(spacing, center_coordinates=True)
# First produce the weighted variance weights
variance_weights = reducer.filter(coordinates, data.velocity_up, weights)[-1]
# And now produce the uncertainty propagation weights
reducer.set_params(uncertainty=True)
block_coords, velocity, uncertainty_weights = reducer.filter(
    coordinates, data.velocity_up, weights

An advantage of using the Green's functions based :class:`verde.Spline` over
:class:`verde.ScipyGridder` is that you can assign weights to the data to incorporate
the data uncertainties or variance into the gridding.
In this example, we'll see how to combine :class:`verde.BlockMean` to decimate the data
and use weights based on the data uncertainty during gridding.
"""
import matplotlib.pyplot as plt
from matplotlib.colors import PowerNorm
import cartopy.crs as ccrs
import pyproj
import numpy as np
import verde as vd

# We'll test this on the California vertical GPS velocity data because it comes with the
# uncertainties
data = vd.datasets.fetch_california_gps()
coordinates = (data.longitude.values, data.latitude.values)

# Use a Mercator projection for our Cartesian gridder
projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())

# Now we can chain a block weighted mean and weighted spline together. We'll use
# uncertainty propagation to calculate the new weights from block mean because our data
# vary smoothly but have different uncertainties.
spacing = 5 / 60  # 5 arc-minutes
chain = vd.Chain(
    [
        ("mean", vd.BlockMean(spacing=spacing * 111e3, uncertainty=True)),
        ("spline", vd.Spline(damping=1e-10)),
    ]
)
print(chain)

can easily weight the data to give each point more or less influence over the results.
This is a good way to not let data points with large uncertainties bias the
interpolation or the data decimation.
"""
# The weights vary a lot so it's better to plot them using a logarithmic color scale
from matplotlib.colors import LogNorm
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import verde as vd

########################################################################################
# We'll use some sample GPS vertical ground velocity which has some variable
# uncertainties associated with each data point. The data are loaded as a
# pandas.DataFrame:
data = vd.datasets.fetch_california_gps()
print(data.head())

########################################################################################
# Let's plot our data using Cartopy to see what the vertical velocities and their
# uncertainties look like. We'll make a function for this so we can reuse it later on.


def plot_data(coordinates, velocity, weights, title_data, title_weights):
    "Make two maps of our data, one with the data and one with the weights/uncertainty"
    fig, axes = plt.subplots(
        1, 2, figsize=(9.5, 7), subplot_kw=dict(projection=ccrs.Mercator())
    )
    crs = ccrs.PlateCarree()
    ax = axes[0]
    ax.set_title(title_data)
    maxabs = vd.maxabs(velocity)

"""
Vector Data
===========

Some datasets have multiple vector components measured for each location, like the East
and West components of wind speed or GPS velocities. For example, let's look at our
sample GPS velocity data from the U.S. West coast.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import pyproj
import verde as vd

data = vd.datasets.fetch_california_gps()

# We'll need to project the geographic coordinates to work with our Cartesian
# classes/functions
projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())
proj_coords = projection(data.longitude.values, data.latitude.values)
# This will be our desired grid spacing in degrees
spacing = 12 / 60

plt.figure(figsize=(6, 8))
ax = plt.axes(projection=ccrs.Mercator())
crs = ccrs.PlateCarree()
tmp = ax.quiver(
    data.longitude.values,
    data.latitude.values,
    data.velocity_east.values,
    data.velocity_north.values,

==================================

Sometimes data has outliers or less reliable points that might skew a blocked mean or
even a median. If the reduction function can take a ``weights`` argument, like
``numpy.average``, you can pass in weights to :class:`verde.BlockReduce` to lower the
influence of the offending data points. However, :class:`verde.BlockReduce` can't
produce weights for the blocked data (for use by a gridder, for example). If you want to
produced blocked weights as well, use :class:`verde.BlockMean`.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import verde as vd

# We'll test this on the California vertical GPS velocity data
data = vd.datasets.fetch_california_gps()

# We'll add some random extreme outliers to the data
outliers = np.random.RandomState(2).randint(0, data.shape[0], size=20)
data.velocity_up[outliers] += 0.08
print("Index of outliers:", outliers)

# Create an array of weights and set the weights for the outliers to a very low value
weights = np.ones_like(data.velocity_up)
weights[outliers] = 1e-5

# Now we can block average the points with and without weights to compare the outputs.
reducer = vd.BlockReduce(reduction=np.average, spacing=30 / 60, center_coordinates=True)
coordinates, no_weights = reducer.filter(
    (data.longitude, data.latitude), data.velocity_up
)
__, with_weights = reducer.filter(