6 - Working with Rasters¶

In this notebook, we will demonstrate how to use the gplately.Raster class to:

  • Create rasters from scattered geographic points
  • Download time-dependent rasters
  • Plot rasters
  • Resize and resample rasters
  • Reconstruct rasters back in time
  • Query rasters using linear interpolation or a spatial tree search
  • Clip rasters by extent
In [1]:
import os, warnings
from pathlib import Path
import cartopy.crs as ccrs
import gplately
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from plate_model_manager import PresentDayRasterManager,PlateModelManager

data_dir = Path("WorkflowData/06-Raster")
data_dir.mkdir(parents=True, exist_ok=True)

Use PlateModelManager to download plate tectonic models, and create PlateReconstruction and PlotTopologies objects.

In [2]:
model_name = "Zahirovic2022"
pmm_model = PlateModelManager().get_model(model_name, data_dir="plate-model-repo")
assert pmm_model is not None, f"Failed to load the {model_name} plate model."
rotation_model = pmm_model.get_rotation_model()
topology_features = pmm_model.get_topologies()
static_polygons = pmm_model.get_static_polygons()

coastlines = pmm_model.get_layer("Coastlines")
continents = pmm_model.get_layer("ContinentalPolygons")
COBs = pmm_model.get_layer("COBs")

model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)
gplot = gplately.PlotTopologies(model, coastlines, continents, COBs)

Alternatively, you can use the gplately.auxiliary.get_gplot() function to get the PlotTopologies and PlateReconstruction objects.

from gplately.auxiliary import get_gplot, get_plate_reconstruction

# use the auxiliary function to create PlotTopologies and PlateReconstruction objects
gplot = get_gplot(model_name) # the PlotTopologies object
model = gplot.plate_reconstruction # the PlateReconstruction object

Create Raster from Scattered Geographic Points¶

In [3]:
with warnings.catch_warnings():
    warnings.filterwarnings("ignore", category=RuntimeWarning)
    import pygmt
    import urllib.request

    relief_filename = data_dir / "earth_relief_10m.nc"
    fallback_url = "https://repo.gplates.org/webdav/gplately/earth_relief_10m.nc"
    if os.path.isfile(relief_filename):
        relief = xr.open_dataarray(relief_filename)
        print(f"Loaded relief data from {relief_filename}")
    else:
        try:
            relief = pygmt.datasets.load_earth_relief(
                resolution="10m", region=[-180, 180, -80, 80]
            )
            print(f"Downloaded relief data from GMT server and saved to {relief_filename}")
            relief.to_netcdf(relief_filename)
        except Exception as err:
            print(f"GMT server download failed ({err}). Trying fallback URL...")
            urllib.request.urlretrieve(fallback_url, relief_filename)
            relief = xr.open_dataarray(relief_filename)
            print(f"Downloaded relief data from fallback URL and saved to {relief_filename}")

    # Randomly subsample the regular grid into "scattered" points.
    # These points will be used to demonstrate the functionality of the Raster.from_points() method later.
    rng = np.random.default_rng(0)
    ny, nx = relief.shape
    n_points = 5000
    iy = rng.integers(0, ny, n_points)
    ix = rng.integers(0, nx, n_points)

    lon = relief.lon.values[ix]
    lat = relief.lat.values[iy]
    val = relief.values[iy, ix]

    # Create a Raster from the points
    raster = gplately.Raster.from_points(
        lon=lon,
        lat=lat,
        values=val,
        spacing="0.5d",  # dataset covers only a small region, so use fine spacing
    )

    fig = pygmt.Figure()
    raster.plot(ax_or_fig=fig, use_gmt=True)
    fig.plot(
        x=lon,
        y=lat,
        style="c0.05c",
        fill="black",
        label="Sample points",
    )
    fig.colorbar(frame='af+l"Earth Relief (m)"')
    fig.legend()
    fig.show()

gmtread [NOTICE]: Remote data courtesy of GMT data server remotedata.generic-mapping-tools.org [http://oceania.generic-mapping-tools.org]

gmtread [NOTICE]: SRTM15 Earth Relief v2.7 at 10x10 arc minutes reduced by Gaussian Cartesian filtering (52.4 km fullwidth) [Tozer et al., 2019].

gmtread [NOTICE]:   -> Download grid file [3.0M]: earth_relief_10m_g.grd

Downloaded relief data from GMT server and saved to WorkflowData/06-Raster/earth_relief_10m.nc
No description has been provided for this image

Download Time-dependent Rasters¶

Let's use PlateModelManager to download netCDF age grids.

There is a unique age grid for each millionth year - let's access the 0 Ma age grid by passing time to the get_raster function.

We can then create a GPlately Raster object from the raster file.

In [4]:
time = 0  # Ma
age_grid_raster = gplately.Raster(
    data=pmm_model.get_raster("AgeGrids", time),
    plate_reconstruction=model,
    extent=(-180, 180, -90, 90),
)
test = gplately.Raster(age_grid_raster)

Plot Raster¶

The age_grid_raster is a Raster object - this object allows us to work with age grids and other rasters. Let's visualise the data with imshow.

In [5]:
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
ax.coastlines() # type: ignore
cpt_file = data_dir / "agegrid.cpt"
if not os.path.isfile(cpt_file):
    import urllib.request

    urllib.request.urlretrieve(
        "https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/tests-dir/unittest/create-age-grids-video/agegrid.cpt",
        cpt_file,
    )
from gplately.plot.gmt_cpt import get_cmap_from_gmt_cpt

age_grid_raster.plot(ax, cmap=get_cmap_from_gmt_cpt(str(cpt_file)))  # Raster.plot() method
ax.set_title("Map Plotted by Cartopy")
plt.show()

/home/runner/micromamba/envs/gplately-all/lib/python3.12/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
No description has been provided for this image

Use the built-in PyGMT plot engine.

In [6]:
with warnings.catch_warnings():
    warnings.filterwarnings("ignore", category=RuntimeWarning)

    from gplately.auxiliary import get_pygmt_basemap_figure

    fig = get_pygmt_basemap_figure(
        projection="N180/10c",
        region="d",
        frame=["xafg30", "yafg30"],
        title="Map Plotted by PyGMT",
    )
    fig.coast(land="darkgreen")
    age_grid_raster.plot(
        fig, cmap=str(cpt_file), nan_transparent=True, use_gmt=True
    )  # Raster.plot() method
    fig.show(crop="+m0.4c")
No description has been provided for this image

Let's plot this netCDF grid along with coastlines, mid-ocean ridges and subduction zones (with teeth).

In [7]:
with warnings.catch_warnings():
    warnings.filterwarnings("ignore", category=UserWarning)

    fig = plt.figure(figsize=(6, 6))
    ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))

    gplot.time = time

    age_grid_raster.plot(ax, cmap=get_cmap_from_gmt_cpt(cpt_file), vmin=0, vmax=200)
    gplot.plot_coastlines(ax, edgecolor="none", facecolor="grey", alpha=0.3)
    gplot.plot_all_topological_sections(
        ax,
        plot_subduction_teeth=True,
        other_kwargs={"color": "grey", "linewidth": 0.8},
        ridge_kwargs={"color": "black", "linewidth": 1.0},
        transform_kwargs={"color": "green", "linewidth": 1.0},
        trench_kwargs={"color": "blue", "linewidth": 1.0},
    )
    ax.set_title(f"{time} Ma")
    plt.show()
No description has been provided for this image

Resize and Resample Raster¶

Let's resize and resample the present-day age grid.

  • resize - provide the number of points in each direction to resize the raster (e.g. 100 cols by 200 rows)
  • resample - provide the grid spacing in each direction (e.g. 0.1 degrees by 0.2 degrees)
In [8]:
age_grid_raster = gplately.Raster(
    data=pmm_model.get_raster("AgeGrids", time),
    extent=(-180, 180, -90, 90),
)

# Set grid size in x and y directions
age_grid_361_181 = age_grid_raster.resize(361, 181)
age_grid_2_degree_space = age_grid_raster.resample(spacingX=2, spacingY=2)

fig = plt.figure(figsize=(16, 8), dpi=100)

ax_1 = fig.add_subplot(221, projection=ccrs.PlateCarree())
ax_1.coastlines()
age_grid_raster.plot(ax=ax_1, cmap=get_cmap_from_gmt_cpt(cpt_file), vmax=200, vmin=0)
ax_1.set_title(f"Original Age Grid {age_grid_raster.shape}")

ax_2 = fig.add_subplot(222, projection=ccrs.PlateCarree())
ax_2.coastlines()
age_grid_361_181.plot(ax=ax_2, cmap=get_cmap_from_gmt_cpt(cpt_file), vmax=200, vmin=0)
ax_2.set_title(f"Age Grid Resized {age_grid_361_181.shape}")

ax_3 = fig.add_subplot(223, projection=ccrs.PlateCarree())
ax_3.coastlines()
age_grid_2_degree_space.plot(
    ax=ax_3, cmap=get_cmap_from_gmt_cpt(cpt_file), vmax=200, vmin=0
)
ax_3.set_title(f"Age Grid Resampled {age_grid_2_degree_space.shape}")

age_grid_raster.resize(91, 46, inplace=True)

ax_4 = fig.add_subplot(224, projection=ccrs.PlateCarree())
ax_4.coastlines()
age_grid_raster.plot(ax=ax_4, cmap=get_cmap_from_gmt_cpt(cpt_file), vmax=200, vmin=0)
ax_4.set_title(f"Age Grid Resized In-Place {age_grid_raster.shape}")

plt.show()
No description has been provided for this image

Download ETOPO1¶

Let's visualise an ETOPO1 relief raster. GPlately also makes it easy to import, using PresentDayRasterManager().get_raster

Alternatively, you can import a local netcdf file using by passing the filename to gplately.Raster or gplately.grids.read_netcdf_grid

In [9]:
with warnings.catch_warnings():
    warnings.filterwarnings("ignore")

    from matplotlib import image

    etopo = gplately.Raster(
        data=image.imread(PresentDayRasterManager().get_raster("ETOPO1_tif"))
    )  # This returns ETOPO1 as a `Raster` object.
    etopo.lats = etopo.lats[::-1]

    print(np.shape(etopo))

(2700, 5400, 3)

etopo is a large (5400 x 2700 pixels) RGB image, with a total of 14,580,000 grid points. We can visualise it using a number of methods, such as Raster.plot.

In [10]:
fig = plt.figure(figsize=(6, 6))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
etopo.plot(ax=ax1, interpolation="none")
ax1.set_title(f"Original ETOPO1")
plt.show()
No description has been provided for this image

We can also resize this RGB raster using the Raster.resize and Raster.resample methods. Here we resample it to a 0.5 by 0.5 degree resolution, then plot it using Raster.plot.

In [11]:
etopo_downscaled = etopo.resample(0.5, 0.5)
print(etopo_downscaled.shape)

fig = plt.figure(figsize=(6, 6))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
ax1.set_title(f"Downsampled ETOPO1")
etopo_downscaled.plot(ax1, interpolation="none")

(361, 721, 3)
Out[11]:
<matplotlib.image.AxesImage at 0x7efc94fb6540>
No description has been provided for this image

Reconstruct Raster¶

The ETOPO1 raster can be reconstructed back in time by assigning a plate_reconstruction to the object. We use the plate reconstruction we defined earlier (model). Note that a Raster object will also accept a plate_reconstruction at initialisation.

After this, we use reconstruct to reconstruct the raster to a given time.

In [12]:
# use the downsampled ETOPO1 so that the reconstruction won't take too long for this demonstration
# assign a plate reconstruction in order to reconstruct the raster
etopo_downscaled.plate_reconstruction = model

white_rgb = (255, 255, 255)  # RGB code for white, to fill gaps in output
etopo_reconstructed = etopo_downscaled.reconstruct(50, threads=4, fill_value=white_rgb)

fig = plt.figure(figsize=(6, 6))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=180))
ax1.set_title(f"Reconstructed ETOPO1 at {etopo_reconstructed.time} Ma")
etopo_reconstructed.plot(ax1)
Out[12]:
<matplotlib.image.AxesImage at 0x7efc94e67020>
No description has been provided for this image

Raster can be reconstructed in-place (inplace=True), and fill_value can be set to any valid matplotlib colour when reconstructing RGB images.

In [13]:
# create a duplicate of the obj
etopo_dup = etopo_downscaled.copy()
etopo_dup.reconstruct(time=75, threads=4, fill_value="darkblue", inplace=True)
etopo_dup.plot(projection=ccrs.Robinson())
plt.gca().set_title(f"Reconstructed ETOPO1 at {etopo_dup.time} Ma")
Out[13]:
Text(0.5, 1.0, 'Reconstructed ETOPO1 at 75.0 Ma')
No description has been provided for this image

By default, Raster.reconstruct uses self.plate_reconstruction.static_polygons to assign plate IDs to grid points. To override this behaviour, pass any collection of pygplates.Feature (e.g. list, pygplates.FeatureCollection, etc) to the partitioning_features argument.

In [14]:
etopo_reconstructed = etopo_downscaled.reconstruct(
    140, partitioning_features=continents, threads=4, fill_value="grey"
)
etopo_reconstructed.plot(projection=ccrs.Orthographic(0, -80))
plt.gca().set_title(f"Reconstructed to {etopo_reconstructed.time} Ma")
Out[14]:
Text(0.5, 1.0, 'Reconstructed to 140.0 Ma')
No description has been provided for this image

Reverse Reconstructions¶

Rasters can also be reverse reconstructed forward in time! Some data will be lost during the reconstruction because that data doesn't exist at 50 Ma.

In [15]:
assert etopo_downscaled.shape[2] == 3, "Must be a RGB image"

etopo_50_ma = etopo_downscaled.reconstruct(50, fill_value="white", threads=4)
etopo_reversed_0_ma = etopo_50_ma.reconstruct(0, fill_value=(0, 0, 0), threads=4)

fig, axs = plt.subplots(
    2,
    2,
    figsize=(10, 5),
    subplot_kw={"projection": ccrs.Mollweide(central_longitude=0)},
)
etopo_downscaled.plot(ax=axs[0][0])
etopo_50_ma.plot(ax=axs[0][1])
etopo_reversed_0_ma.plot(ax=axs[1][0])

use_spatial_tree = True  # set this to False to use Linear Interpolation to fill the gaps. (much slower!!)
if use_spatial_tree:
    etopo_reversed_0_ma.fill_gaps(invalid_value=(0, 0, 0), use_spatial_tree=True).plot(
        ax=axs[1][1]
    )
else:
    etopo_reversed_0_ma.fill_gaps(invalid_value=(0, 0, 0), method="linear").plot(
        ax=axs[1][1]
    )

axs[0][0].set_title("Original Present Day")
axs[0][1].set_title("Reconstructed to 50 Ma")
axs[1][0].set_title("Reverse from 50 Ma to Present Day")
if use_spatial_tree:
    axs[1][1].set_title("Gaps Filled By Nearest Interpolation")
else:
    axs[1][1].set_title("Gaps Filled By Linear Interpolation")
plt.show()
No description has been provided for this image

Reconstructing NetCDF Raster¶

Similar to above, we can also reconstruct numeric netcdf grids! GPlately also includes ETOPO1 as a netcdf for download.

In [16]:
# download ETOPO1 netCDF
etopo_nc = gplately.Raster(data=PresentDayRasterManager().get_raster("ETOPO1_grd"))
etopo_nc._data = etopo_nc._data.astype(float)

# resample to a more manageable size
etopo_nc.resample(0.5, 0.5, inplace=True)
print(etopo_nc.shape)

# Assign plate reconstruction
etopo_nc.plate_reconstruction = model

# Reconstruct raster to 50 Ma
etopo_nc_reconstructed = etopo_nc.reconstruct(50, threads=4)

# plot
fig = plt.figure(figsize=(10, 10))
ax_1 = fig.add_subplot(121, projection=ccrs.Mollweide(central_longitude=20))
im = etopo_nc_reconstructed.plot(
    ax_1, cmap=gplately.get_topo_cmap(), vmin=-10927, vmax=8726
)
fig.colorbar(im, ax=ax_1, pad=0.05, shrink=0.7, orientation="horizontal")
ax_1.set_title(f"ETOPO1 NetCDF at {etopo_nc_reconstructed.time} Ma")

ax_2 = fig.add_subplot(122, projection=ccrs.Mollweide(central_longitude=20))
im = etopo_nc_reconstructed.fill_gaps(use_gmt=True).plot(
    ax_2, cmap=gplately.get_topo_cmap(), vmin=-10927, vmax=8726
)
fig.colorbar(im, ax=ax_2, pad=0.05, shrink=0.7, orientation="horizontal")
ax_2.set_title(f"Gaps Filled by Nearest using pygmt.grdfill()")

plt.show()

(361, 721)
No description has been provided for this image
In [17]:
# Save the reconstructed ETOPO grid to a netCDF file
etopo_nc_reconstructed.save_to_netcdf4(
    data_dir / f"reconstructed_etopo_{etopo_nc_reconstructed.time}.nc",
)

Raster Query by Linear Interpolation¶

The code cell below demonstrates how to use Raster.query() with linear interpolation to query data from a raster.

In [18]:
fig = plt.figure(figsize=(12, 3), dpi=100)

ax_1 = fig.add_subplot(131, projection=ccrs.PlateCarree(central_longitude=0))
etopo_downscaled.plot(ax=ax_1)

extent = (-100, 100, -50, 50)
ax_2 = fig.add_subplot(132, projection=ccrs.PlateCarree(central_longitude=0))
xx, yy = np.meshgrid(
    np.linspace(extent[0], extent[1], 100), np.linspace(extent[2], extent[3], 50)
)
shape = xx.shape
xx = xx.flatten()
yy = yy.flatten()
values = etopo_downscaled.query(lons=xx, lats=yy, interpolation_method="linear")
ax_2.set_extent(extent, crs=ccrs.PlateCarree())
ax_2.set_facecolor("white")
ax_2.scatter(
    xx,
    yy,
    c=values / 255.0,
    marker="o",
    s=3.5,
    transform=ccrs.PlateCarree(),
)

ax_3 = fig.add_subplot(133, projection=ccrs.PlateCarree(central_longitude=0))
ax_3.set_global()
ax_3.imshow(
    np.flipud(np.reshape(values, (shape[0], shape[1], 3))),
    transform=ccrs.PlateCarree(),
    extent=extent,
)

ax_1.set_title(f"Original Image")
ax_2.set_title(f"Values Plotted by scatter()")
ax_3.set_title(f"Values Plotted by imshow()")
fig.suptitle("Raster Query by Linear Interpolation")
fig.tight_layout()
plt.show()
No description has been provided for this image

This example demonstrates how to query nearest values from a Raster. We'll find which subduction zones form continental arcs by using a continental mask raster, and querying that raster with points projected 250 km from the trench (in the direction of the subducting plate) to determine if they are inside a continent.

In [19]:
# get the continental mask raster
continental_mask_file = os.path.join("NotebookFiles", "continental_grid_0.nc")
if not os.path.isfile(continental_mask_file):
    import urllib.request

    urllib.request.urlretrieve(
        "https://github.com/GPlates/gplately/raw/refs/heads/master/Notebooks/NotebookFiles/continental_grid_0.nc",
        continental_mask_file,
    )

continental_raster = gplately.Raster(continental_mask_file, model)

# tessellate trenches and extract the subduction polarity angle and the lat-lon coordinates
trench_data = model.tessellate_subduction_zones(time)
trench_normal_azimuthal_angle = trench_data[:, 7]
trench_pt_lon = trench_data[:, 0]
trench_pt_lat = trench_data[:, 1]

# calculate 250 km arc distance
arc_distance = 250 / (gplately.tools.geocentric_radius(trench_pt_lat) / 1e3)

# Lat and lon coordinates of all trench points after being projected out 250 km in the direction of subduction.
dlon = arc_distance * np.sin(np.radians(trench_normal_azimuthal_angle))
dlat = arc_distance * np.cos(np.radians(trench_normal_azimuthal_angle))
ilon = trench_pt_lon + np.degrees(dlon)
ilat = trench_pt_lat + np.degrees(dlat)

Now use these projected trench points to query the nearest cells in the continental raster using Raster.query.

In [20]:
# Query the raster with the projected trench points
sampled_points = continental_raster.query(
    lons=ilon, lats=ilat, interpolation_method="nearest"
)

# The cells of land in the continental raster are 1. Now find all the inland points.
in_raster_indices = sampled_points > 0

# Get the lat-lon coordinates of the in_raster points
lat_in = ilat[in_raster_indices]
lon_in = ilon[in_raster_indices]

Plot the inland points along with the raster, trenches, and coastlines.

In [21]:
fig = plt.figure(figsize=(8, 6), dpi=100)
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
gplot.time = time
gplot.plot_grid_from_netCDF(
    ax1, continental_mask_file, cmap="twilight", alpha=0.5, vmin=0, vmax=200
)
gplot.plot_coastlines(ax1, edgecolor="k", facecolor="1", alpha=0.1)
gplot.plot_trenches(ax1, color="r", zorder=5)  # Plot the original trench points in Red

# Plot the projected inland points in Blue
ax1.plot(
    lon_in,
    lat_in,
    linestyle="none",
    marker="o",
    markersize=0.25,
    markerfacecolor="blue",
    markeredgecolor="blue",
    transform=ccrs.PlateCarree(),
)
plt.title(f"{time} Ma")

# Custom legend
from matplotlib.lines import Line2D

handles = [
    Line2D([0], [0], color="red"),
    Line2D([0], [0], color="blue"),
]
labels = [
    "Trenches",
    "Arc segments in \ncontinental grids",
]

plt.legend(handles, labels, loc="lower left", bbox_to_anchor=(0.0, -0.05))
plt.show()
No description has been provided for this image

Query Raster by Spatial Tree¶

We will reuse the inland sample points from previous code cells to demonstrate how to query a raster within a region of interest. This is slower than interpolation. However, this is useful if you need to find the nearest valid data points within a region of interest.

In [22]:
age_grid_raster = gplately.Raster(data=pmm_model.get_raster("AgeGrids", 0))

# plot the age grid raster and black sample points
fig = plt.figure(figsize=(10, 10), dpi=100)
ax_1 = fig.add_subplot(121, projection=ccrs.Mollweide(central_longitude=0))
age_grid_raster.plot(ax=ax_1, cmap=get_cmap_from_gmt_cpt(cpt_file), vmax=200, vmin=0)
ax_1.plot(
    lon_in,
    lat_in,
    linestyle="none",
    marker="o",
    markersize=0.25,
    markerfacecolor="black",
    markeredgecolor="black",
    transform=ccrs.PlateCarree(),
)
ax_1.set_title("Age Grid and Black Sample Points")

# plot the data being retrieved by raster query
ax_2 = fig.add_subplot(122, projection=ccrs.Mollweide(central_longitude=0))
ax_2.set_global()
xx = lon_in
yy = lat_in

# You may set the `region_of_interest` to a smaller value, such as 100, to see what will happen.
# You will see fewer points being plotted because there are no valid data within 100 km for some sample points.
roi = 500  # km
values = age_grid_raster.query(lons=xx, lats=yy, region_of_interest=roi)

ax_2.scatter(
    xx,
    yy,
    c=values,
    marker="o",
    s=0.25,
    transform=ccrs.PlateCarree(),
    cmap=get_cmap_from_gmt_cpt(cpt_file),
    vmax=200,
    vmin=0,
)

gl = ax_2.gridlines(
    crs=ccrs.PlateCarree(),
    draw_labels=False,
    linewidth=1,
    color="gray",
    alpha=0.5,
    linestyle="--",
)
ax_2.set_title(f"Sampled Age Values within {roi} KM")
fig.tight_layout()
plt.show()
No description has been provided for this image

Clip Raster by Extent¶

In [23]:
fig, axs = plt.subplots(
    1,
    2,
    figsize=(10, 8),
    gridspec_kw={"width_ratios": [2, 1]},
    subplot_kw={"projection": ccrs.PlateCarree()},
)
# fig.tight_layout()
ax_1 = axs[0]
ax_2 = axs[1]

# plot the original age grid raster
age_grid_raster.plot(ax=ax_1, cmap=get_cmap_from_gmt_cpt(cpt_file), vmax=200, vmin=0)
ax_1.set_title("Original Age Grid")
gl = ax_1.gridlines(
    crs=ccrs.PlateCarree(),
    draw_labels=True,
    linewidth=1,
    color="gray",
    alpha=0.5,
    linestyle="--",
)
gl.right_labels = False
gl.top_labels = False

# clip the raster with a new extent
clipped_raster = age_grid_raster.clip_by_extent((-50, 50, -80, 40))

# plot the clipped raster
clipped_raster.plot(
    ax=ax_2,
    transform=ccrs.PlateCarree(),
    cmap=get_cmap_from_gmt_cpt(cpt_file),
    vmax=200,
    vmin=0,
)
ax_2.set_title("Clipped Age Grid")
gl = ax_2.gridlines(
    crs=ccrs.PlateCarree(),
    draw_labels=True,
    linewidth=1,
    color="gray",
    alpha=0.5,
    linestyle="--",
)
gl.left_labels = False
gl.top_labels = False
plt.show()
No description has been provided for this image