How to use the verde.grid_coordinates function in verde

# Now we can use the bi-harmonic spline method [Sandwell1987]_ to fit this data. First,
# we create a new :class:`verde.Spline`:

spline = vd.Spline()
# Printing a gridder shows the class and all of it's configuration options.
print(spline)

########################################################################################
# Before we can use the spline, we need to fit it to our synthetic data. After that, we
# can use the spline to predict values anywhere:

spline.fit((data.easting, data.northing), data.scalars)

# Generate coordinates for a regular grid with 100 m grid spacing (assuming coordinates
# are in meters).
grid_coords = vd.grid_coordinates(region=(0, 5000, -5000, 0), spacing=100)
gridded_scalars = spline.predict(grid_coords)

plt.figure()
plt.pcolormesh(grid_coords[0], grid_coords[1], gridded_scalars, cmap="RdBu_r")
plt.colorbar()
plt.show()

########################################################################################
# We can compare our predictions with the true values for the checkerboard function
# using the :meth:`~verde.Spline.score` method to calculate the `R² coefficient of
# determination `__.

true_values = vd.datasets.CheckerBoard().predict(grid_coords)
print(spline.score(grid_coords, true_values))

########################################################################################

import boule as bl
import harmonica as hm


# Use the WGS84 ellipsoid to obtain the mean Earth radius which we'll use to
# reference the tesseroid
ellipsoid = bl.WGS84

# Define tesseroid with top surface at the mean Earth radius, thickness of 10km
# (bottom = top - thickness) and density of 2670kg/m^3
tesseroid = [-70, -50, -40, -20, ellipsoid.mean_radius - 10e3, ellipsoid.mean_radius]
density = 2670

# Define computation points on a regular grid at 100km above the mean Earth
# radius
coordinates = vd.grid_coordinates(
    region=[-80, -40, -50, -10],
    shape=(80, 80),
    extra_coords=100e3 + ellipsoid.mean_radius,
)

# Compute the radial component of the acceleration
gravity = hm.tesseroid_gravity(coordinates, tesseroid, density, field="g_z")
print(gravity)

# Plot the gravitational field
fig = plt.figure(figsize=(8, 9))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-60))
img = ax.pcolormesh(*coordinates[:2], gravity, transform=ccrs.PlateCarree())
plt.colorbar(img, ax=ax, pad=0, aspect=50, orientation="horizontal", label="mGal")
ax.coastlines()
ax.set_title("Downward component of gravitational acceleration")

Normal gravity is defined as the magnitude of the gradient of the gravity
potential (gravitational + centrifugal) of a reference ellipsoid. Function
:func:`harmonica.normal_gravity` implements a closed form solution
[LiGotze2001]_ to calculate normal gravity at any latitude and height (it's
symmetric in the longitudinal direction). The ellipsoid in the calculations
used can be changed using :func:`harmonica.set_ellipsoid`.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import harmonica as hm
import verde as vd

# Create a global computation grid with 1 degree spacing
region = [0, 360, -90, 90]
longitude, latitude = vd.grid_coordinates(region, spacing=1)

# Compute normal gravity for the WGS84 ellipsoid (the default) on its surface
gamma = hm.normal_gravity(latitude=latitude, height=0)

# Make a plot of the normal gravity using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic())
ax.set_title("Normal gravity of the WGS84 ellipsoid")
ax.coastlines()
pc = ax.pcolormesh(longitude, latitude, gamma, transform=ccrs.PlateCarree())
plt.colorbar(
    pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.5
)
plt.tight_layout()
plt.show()

plt.xlabel("Easting")
plt.ylabel("Northing")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, 1.18))
plt.show()

######################################################################################
# Extra Coordinates
# -----------------
#
# In some cases, you might need an additional coordinate such as a height or a time
# that is associated with your coordinate grid. The ``extra_coords`` parameter
# in :func:`verde.grid_coordinates` creates an extra coordinate array that is the same
# shape as the coordinate grid, but contains a constant value. For example, let's
# add a constant height of 1000 units and time of 1 to our coordinate grid.

easting, northing, height, time = vd.grid_coordinates(
    region=region, spacing=spacing, extra_coords=[1000, 1]
)

print(easting.shape, northing.shape, height.shape, time.shape)

########################################################################################
# And we can print the height array to verify that it is correct

print(height)

########################################################################################
# And we can print the time array as well

print(time)

import numpy as np
import matplotlib.pyplot as plt
import verde as vd
from matplotlib.patches import Rectangle

########################################################################################
# First let's create a region that is 1000 units west-east and 1000 units south-north,
# and we will set an initial spacing to 100 units.

spacing = 100
west, east, south, north = 0, 1000, 0, 1000
region = (west, east, south, north)

# create the grid coordinates
easting, northing = vd.grid_coordinates(region=region, spacing=spacing)

########################################################################################
# We can check the dimensions of the grid coordinates. The region is 1000 units and the
# spacing is 100 units, so the shape of the segments is 10x10. However, the number of
# grid nodes in this case is one more than the number of segments. So our grid
# coordinates have a shape of 11x11.
print(easting.shape, northing.shape)

########################################################################################
# Let's define two functions to visualize the region bounds and grid points


def plot_region(ax, region):
    "Plot the region as a solid line."
    west, east, south, north = region
    ax.add_patch(

########################################################################################
# Let's plot our grid side-by-side with one generated by the default spline:

grid_default = spline_default.grid(
    region=region,
    spacing=spacing,
    projection=projection,
    dims=["latitude", "longitude"],
    data_names=["temperature"],
)

mask = vd.distance_mask(
    (data.longitude, data.latitude),
    maxdist=3 * spacing * 111e3,
    coordinates=vd.grid_coordinates(region, spacing=spacing),
    projection=projection,
)

grid = grid.where(mask)
grid_default = grid_default.where(mask)

plt.figure(figsize=(14, 8))
for i, title, grd in zip(range(2), ["Defaults", "Tuned"], [grid_default, grid]):
    ax = plt.subplot(1, 2, i + 1, projection=ccrs.Mercator())
    ax.set_title(title)
    pc = grd.temperature.plot.pcolormesh(
        ax=ax,
        cmap="plasma",
        transform=ccrs.PlateCarree(),
        vmin=data.air_temperature_c.min(),
        vmax=data.air_temperature_c.max(),

"""
Gravity forward modeling with prisms
====================================

Meh.
"""
import matplotlib.pyplot as plt
import verde as vd
import harmonica as hm

prism = (-500, 500, -1000, 1000, 0, 500)
coordinates = vd.grid_coordinates(
    (-2e3, 2e3, -2e3, 2e3), spacing=10, extra_coords=[-20]
)

gravity = hm.prism_gravity(coordinates, prism, density=2670, field="gz")
print(coordinates, gravity)

plt.figure()
plt.pcolormesh(*coordinates[:2], gravity)
plt.colorbar()
plt.axis("scaled")
plt.show()

########################################################################################
# Adjusting region boundaries when creating the grid coordinates
# --------------------------------------------------------------
#
# Now let's change our spacing to 300 units. Because the range of the west-east and
# south-north boundaries are not multiples of 300, we must choose to change either:
#
# - the boundaries of the region in order to fit the spacing, or
# - the spacing in order to fit the region boundaries.
#
# We could tell :func:`verde.grid_coordinates` to adjust the region boundaries by
# passing ``adjust="region"``.

spacing = 300
region_easting, region_northing = vd.grid_coordinates(
    region=region, spacing=spacing, adjust="region"
)
print(region_easting.shape, region_northing.shape)

########################################################################################
# With the spacing set at 300 units and a 4 by 4 grid of regular dimensions,
# :func:`verde.grid_coordinates` calculates the spatial location of each
# grid point and adjusts the region so that the maximum northing and maximum
# easting values are divisible by the spacing. In this example, both the  easting and
# northing have 3 segments (4 nodes) that are each 300 units long, meaning the easting
# and northing span from 0 to 900. Both dimensions are divisible
# by 300.

print(region_easting)
print(region_northing)