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

# Plot the data uncertainties
ax = axes[0]
ax.set_title("Data uncertainty")
# Plot the uncertainties in mm/yr and using a power law for the color scale to highlight
# the smaller values
pc = ax.scatter(
    *coordinates,
    c=data.std_up * 1000,
    s=20,
    cmap="magma",
    transform=crs,
    norm=PowerNorm(gamma=1 / 2)
)
cb = plt.colorbar(pc, ax=ax, orientation="horizontal", pad=0.05)
cb.set_label("uncertainty [mm/yr]")
vd.datasets.setup_california_gps_map(ax, region=region)
# Plot the gridded velocities
ax = axes[1]
ax.set_title("Weighted spline interpolated velocity")
maxabs = vd.maxabs(data.velocity_up) * 1000
pc = (grid.velocity * 1000).plot.pcolormesh(
    ax=ax,
    cmap="seismic",
    vmin=-maxabs,
    vmax=maxabs,
    transform=crs,
    add_colorbar=False,
    add_labels=False,
)
cb = plt.colorbar(pc, ax=ax, orientation="horizontal", pad=0.05)
cb.set_label("vertical velocity [mm/yr]")
ax.scatter(*coordinates, c="black", s=0.5, alpha=0.1, transform=crs)

cb = plt.colorbar(tmp, ax=ax, orientation="horizontal", pad=0.05)
    cb.set_label("meters/year")
    # Plot the original data
    ax.quiver(
        data.longitude.values,
        data.latitude.values,
        data.velocity_east.values,
        data.velocity_north.values,
        scale=0.2,
        transform=crs,
        color="k",
        width=0.001,
        label="Original data",
    )
    # Setup the map ticks
    vd.datasets.setup_california_gps_map(
        ax, land=None, ocean=None, region=vd.get_region((data.longitude, data.latitude))
    )
    ax.coastlines(color="white")
axes[0].legend(loc="lower left")
plt.tight_layout()
plt.show()

ax = axes[1]
ax.set_title("Weighted spline interpolated velocity")
maxabs = vd.maxabs(data.velocity_up) * 1000
pc = (grid.velocity * 1000).plot.pcolormesh(
    ax=ax,
    cmap="seismic",
    vmin=-maxabs,
    vmax=maxabs,
    transform=crs,
    add_colorbar=False,
    add_labels=False,
)
cb = plt.colorbar(pc, ax=ax, orientation="horizontal", pad=0.05)
cb.set_label("vertical velocity [mm/yr]")
ax.scatter(*coordinates, c="black", s=0.5, alpha=0.1, transform=crs)
vd.datasets.setup_california_gps_map(ax, region=region)
ax.coastlines()
plt.tight_layout()
plt.show()

label="Outliers",
    )
    # Plot the block means and saturate the colorbar a bit to better show the
    # differences in the data.
    pc = ax.scatter(
        *coordinates,
        c=velocity,
        s=70,
        transform=crs,
        cmap="seismic",
        vmin=-maxabs / 3,
        vmax=maxabs / 3
    )
    cb = plt.colorbar(pc, ax=ax, orientation="horizontal", pad=0.05)
    cb.set_label("vertical velocity [m/yr]")
    vd.datasets.setup_california_gps_map(ax)
ax.legend(loc="lower left")
plt.tight_layout()
plt.show()

)

plt.figure(figsize=(6, 8))
ax = plt.axes(projection=ccrs.Mercator())
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.2, 0.15, 0.05, label="0.05 m/yr", coordinates="figure")
ax.set_title("Gridded velocities")
vd.datasets.setup_california_gps_map(ax)
plt.tight_layout()
plt.show()

)
    crs = ccrs.PlateCarree()
    ax = axes[0]
    ax.set_title(title_data)
    maxabs = vd.maxabs(velocity)
    pc = ax.scatter(
        *coordinates,
        c=velocity,
        s=30,
        cmap="seismic",
        vmin=-maxabs,
        vmax=maxabs,
        transform=crs,
    )
    plt.colorbar(pc, ax=ax, orientation="horizontal", pad=0.05).set_label("m/yr")
    vd.datasets.setup_california_gps_map(ax)
    ax = axes[1]
    ax.set_title(title_weights)
    pc = ax.scatter(
        *coordinates, c=weights, s=30, cmap="magma", transform=crs, norm=LogNorm()
    )
    plt.colorbar(pc, ax=ax, orientation="horizontal", pad=0.05)
    vd.datasets.setup_california_gps_map(ax)
    plt.tight_layout()
    plt.show()

crs = ccrs.PlateCarree()
fig, axes = plt.subplots(
    1, 2, figsize=(8, 4), subplot_kw=dict(projection=ccrs.Mercator())
)
# Plot the horizontal velocity vectors
ax = axes[0]
ax.set_title("GPS horizontal velocities")
ax.quiver(
    data.longitude.values,
    data.latitude.values,
    data.velocity_east.values,
    data.velocity_north.values,
    scale=0.3,
    transform=crs,
)
vd.datasets.setup_california_gps_map(ax)
# Plot the vertical velocity
ax = axes[1]
ax.set_title("Vertical velocity")
maxabs = vd.maxabs(data.velocity_up)
tmp = ax.scatter(
    data.longitude,
    data.latitude,
    c=data.velocity_up,
    s=10,
    vmin=-maxabs / 3,
    vmax=maxabs / 3,
    cmap="seismic",
    transform=crs,
)
plt.colorbar(tmp, ax=ax).set_label("meters/year")
vd.datasets.setup_california_gps_map(ax)

# Plot the original data uncertainties
ax = axes[2]
ax.set_title("Data uncertainty")
# Use a power law for the color scale because there is a lot of variability
# Convert m/year to mm/year to have smaller values on the color bar
pc = ax.scatter(
    *coordinates,
    c=data.std_up * 1000,
    s=10,
    transform=crs,
    alpha=1,
    cmap="magma",
    norm=PowerNorm(gamma=1 / 2)
)
plt.colorbar(pc, ax=ax, orientation="horizontal", pad=0.05).set_label("mm/yr")
vd.datasets.setup_california_gps_map(ax)
plt.tight_layout()
plt.show()

components = [grid.east_component, grid.north_component]
for ax, component, title in zip(axes, components, titles):
    ax.set_title(title)
    maxabs = vd.maxabs(component)
    tmp = component.plot.pcolormesh(
        ax=ax,
        vmin=-maxabs,
        vmax=maxabs,
        cmap="bwr",
        transform=crs,
        add_colorbar=False,
        add_labels=False,
    )
    cb = plt.colorbar(tmp, ax=ax, orientation="horizontal", pad=0.05)
    cb.set_label("meters/year")
    vd.datasets.setup_california_gps_map(ax, land=None, ocean=None)
    ax.coastlines(color="white")
plt.tight_layout()
plt.show()

########################################################################################
# Gridding
# --------
#
# You can use :class:`verde.Vector` to create multi-component gridders out of
# :class:`verde.Spline` the same way as we did for trends. In this case, each component
# is treated separately.
#
# We can start by splitting the data into training and testing sets (see
# :ref:`model_selection`). Notice that :func:`verde.train_test_split` work for
# multicomponent data automatically.

# Plot the original data uncertainties
ax = axes[2]
ax.set_title("Data uncertainty")
# Use a power law for the color scale because there is a lot of variability
# Convert m/year to mm/year to have smaller values on the color bar
pc = ax.scatter(
    *coordinates,
    c=data.std_up * 1000,
    s=10,
    transform=crs,
    alpha=1,
    cmap="magma",
    norm=PowerNorm(gamma=1 / 2)
)
plt.colorbar(pc, ax=ax, orientation="horizontal", pad=0.05).set_label("mm/yr")
vd.datasets.setup_california_gps_map(ax)
plt.tight_layout()
plt.show()