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

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
from importlib.resources import files
from gplately.plot import get_age_grid_cmap

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

Matplotlib is building the font cache; this may take a moment.

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

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

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 [ERROR]: Remote download is currently deactivated

gmtread [ERROR]: Unable to obtain remote file @earth_relief_10m_g

gmtread [ERROR]: Remote download is currently deactivated

gmtread [ERROR]: Unable to obtain remote file @earth_relief_10m_g

gmtread [ERROR]: Remote download is currently deactivated

gmtread [ERROR]: Unable to obtain remote file @earth_relief_10m_g

gmtread [ERROR]: File @earth_relief_10m_g not found

[Session pygmt-session (3)]: Error returned from GMT API: GMT_FILE_NOT_FOUND (16)

GMT server download failed (Module 'read' failed with status code 16:
gmtread [ERROR]: Remote download is currently deactivated
gmtread [ERROR]: Unable to obtain remote file @earth_relief_10m_g
gmtread [ERROR]: Remote download is currently deactivated
gmtread [ERROR]: Unable to obtain remote file @earth_relief_10m_g
gmtread [ERROR]: Remote download is currently deactivated
gmtread [ERROR]: Unable to obtain remote file @earth_relief_10m_g
gmtread [ERROR]: File @earth_relief_10m_g not found). Trying fallback URL...

Downloaded relief data from fallback URL and saved to WorkflowData/06-Raster/earth_relief_10m.nc

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.

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.

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
ax.coastlines() # type: ignore
age_grid_raster.plot(ax, cmap=get_age_grid_cmap()) 
ax.set_title("Map Plotted by Cartopy")
plt.show()

/Users/runner/micromamba/envs/gplately-notebook-env/lib/python3.14/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip

Use the built-in PyGMT plot engine.

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(files("gplately").joinpath("data", "agegrid.cpt")), nan_transparent=True, use_gmt=True
    )  
    fig.show(crop="+m0.4c")

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

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_age_grid_cmap(), 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()

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)

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_age_grid_cmap(), 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_age_grid_cmap(), 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_age_grid_cmap(), 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_age_grid_cmap(), vmax=200, vmin=0)
ax_4.set_title(f"Age Grid Resized In-Place {age_grid_raster.shape}")

plt.show()

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

with warnings.catch_warnings():
    warnings.filterwarnings("ignore")

    from matplotlib import image

    etopo = gplately.Raster(
        data=image.imread(PresentDayRasterManager().get_raster("ETOPO1_tif")),
        origin="upper",
    )

    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.

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()

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.

etopo_downscaled = etopo.resample(0.5, 0.5)
print(etopo_downscaled.shape)
print(etopo_downscaled.dtype)

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)
uint8

<matplotlib.image.AxesImage at 0x159bfb250>

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.

# 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)

<matplotlib.image.AxesImage at 0x150f93390>

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

# 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")

Text(0.5, 1.0, 'Reconstructed ETOPO1 at 75.0 Ma')

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.

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")

Text(0.5, 1.0, 'Reconstructed to 140.0 Ma')

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.

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()

Reconstructing NetCDF Raster

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

# 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)

# 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.

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()

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.

# 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.

# 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.

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()

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.

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_age_grid_cmap(), 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_age_grid_cmap(),
    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()

Clip Raster

from gplately.utils.feature_filter import (
    FeatureNameFilter,
    PolygonAreaFilter,
    filter_feature_collection,
)

# find the polygon feature with area greater than 7 million and named "Australia",
# which is the mainland of Australia, and use the Australia mainland polygon to clip a raster.
features = filter_feature_collection(
    gplately.load_feature_collection(pmm_model.get_coastlines()),  # type: ignore
    [
        FeatureNameFilter(
            ["Australia"],
            exact_match=True,
            case_sensitive=True,
        ),
        PolygonAreaFilter(7e6, reverse=False),
    ],
)
assert (
    len(features) == 1
), f"Expected to find exactly one feature for the mainland of Australia, but found {len(features)} features."

australia_mainland_geometry = features[0].get_geometry()
# now the austrlia is in a global map
raster_australia = etopo_downscaled.clip_by_polygons(
    [australia_mainland_geometry], fill_value="lightblue"
)

lat_lon = australia_mainland_geometry.to_lat_lon_array()
min_lat = np.min(lat_lon[:, 0])
max_lat = np.max(lat_lon[:, 0])
min_lon = np.min(lat_lon[:, 1])
max_lon = np.max(lat_lon[:, 1])
australia_extent = (
    int(np.floor(min_lon)) - 1,
    int(np.ceil(max_lon)) + 1,
    int(np.floor(min_lat)) - 1,
    int(np.ceil(max_lat)) + 1,
)

clipped_australia = raster_australia.clip_by_extent(australia_extent)

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

# plot the original raster
ax_1.set_global()
etopo_downscaled.plot(ax=ax_1)
gl = ax_1.gridlines(
    draw_labels=False,
    linewidth=1,
    color="gray",
    alpha=0.5,
    linestyle="--",
)
gl.right_labels = False
gl.top_labels = False

# plot the clipped raster
ax_2.set_global()
raster_australia.plot(ax=ax_2)
gl = ax_2.gridlines(
    draw_labels=False,
    linewidth=1,
    color="gray",
    alpha=0.5,
    linestyle="--",
)

# plot the clipped raster with tighter extent
ax_3.set_extent(australia_extent)
clipped_australia.plot(ax=ax_3)
gl = ax_3.gridlines(
    draw_labels=False,
    linewidth=1,
    color="gray",
    alpha=0.5,
    linestyle="--",
)
gl.right_labels = True
gl.bottom_labels = True

plt.show()

2026-09-14 01:37:31 - gplately - WARNING - The new data shape (66, 87, 3) is different from the old shape (361, 721, 3)