How do I control the size of margins around a cartopy map projection?

I'm trying to plot a bunch of data on a map of the sky in various projections, using matplotlib + cartopy, and the margins around the maps are always too large and none of the controls I can find seem to help. Example (annotations added after rendering):

I would like to make the outer margin of the entire image, including the colorbar, be something like 2.5mm, and the gap between the colorbar and the image be something like 5mm (these numbers will need to be tweaked of course), and then the map should fill the rest of the available space.

Note that I may need to turn this into a 2-subplot figure, two maps sharing a colorbar, each with a label, and possibly also add meridians and parallels with labels, so a solution that works regardless of how much 'furniture' each Axes has is strongly preferred.

Part of the problem seems to be that each map projection has its own desired aspect ratio, and if that doesn't agree with the aspect ratio of the figure then spacing will be added to preserve said aspect ratio, but since the desired aspect ratio is not documented anywhere and the width of the colorbar is unpredictable, knowing that doesn't actually help me any. Also, this really is only part of the problem; if I hold the overall figure height constant and vary the width over a range of values, the figure has the least amount of unwanted white space when the figure's aspect ratio is just so, but it still has unwanted white space.

Here's one version of the code I have now. Please note how each projection gets rendered with different margins.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.layout_engine import ConstrainedLayoutEngine


def main():
    cel_sphere = ccrs.Globe(datum=None, ellipse=None,
                            semimajor_axis=180/np.pi,
                            semiminor_axis=180/np.pi)
    sky_plate = ccrs.PlateCarree(globe=cel_sphere)

    ra, dec = np.mgrid[-179.5:180:1, -89.5:90:1]
    fake_data = generate_perlin_noise_2d(ra.shape, (1, 1))

    for label, proj in [("ee", ccrs.EqualEarth),
                        ("mw", ccrs.Mollweide),
                        ("lc", ccrs.LambertCylindrical)]:
        try:
            fig, ax = plt.subplots(
                figsize=(20, 10),
                layout=ConstrainedLayoutEngine(
                    h_pad=0, w_pad=0, hspace=0, wspace=0
                ),
                subplot_kw={
                    "xlim": (-180, 180),
                    "ylim": (-90, 90),
                    "projection": proj(globe=cel_sphere)
                },
            )
            ctr = ax.contourf(ra, dec, fake_data,
                              transform=sky_plate,
                              cmap="Greys")
            fig.colorbar(ctr, shrink=0.5, pad=0.02)
            fig.savefig(f"layout_test_{label}.png")
        finally:
            plt.close(fig)

# stolen from https://pvigier.github.io/2018/06/13/perlin-noise-numpy.html
def generate_perlin_noise_2d(shape, res):
    def f(t):
        return 6*t**5 - 15*t**4 + 10*t**3

    delta = (res[0] / shape[0], res[1] / shape[1])
    d = (shape[0] // res[0], shape[1] // res[1])
    grid = np.mgrid[0:res[0]:delta[0],0:res[1]:delta[1]].transpose(1, 2, 0) % 1
    # Gradients
    angles = 2*np.pi*np.random.rand(res[0]+1, res[1]+1)
    gradients = np.dstack((np.cos(angles), np.sin(angles)))
    g00 = gradients[0:-1,0:-1].repeat(d[0], 0).repeat(d[1], 1)
    g10 = gradients[1:,0:-1].repeat(d[0], 0).repeat(d[1], 1)
    g01 = gradients[0:-1,1:].repeat(d[0], 0).repeat(d[1], 1)
    g11 = gradients[1:,1:].repeat(d[0], 0).repeat(d[1], 1)
    # Ramps
    n00 = np.sum(grid * g00, 2)
    n10 = np.sum(np.dstack((grid[:,:,0]-1, grid[:,:,1])) * g10, 2)
    n01 = np.sum(np.dstack((grid[:,:,0], grid[:,:,1]-1)) * g01, 2)
    n11 = np.sum(np.dstack((grid[:,:,0]-1, grid[:,:,1]-1)) * g11, 2)
    # Interpolation
    t = f(grid)
    n0 = n00*(1-t[:,:,0]) + t[:,:,0]*n10
    n1 = n01*(1-t[:,:,0]) + t[:,:,0]*n11
    return np.sqrt(2)*((1-t[:,:,1])*n0 + t[:,:,1]*n1)

main()

And here's a version that renders two subplots with all possible labels, demonstrating that the unwanted space is not just because of leaving space for furniture.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np

from matplotlib import colors, cm
from matplotlib.layout_engine import ConstrainedLayoutEngine


def main():
    cel_sphere = ccrs.Globe(datum=None, ellipse=None,
                            semimajor_axis=180/np.pi,
                            semiminor_axis=180/np.pi)
    sky_plate = ccrs.PlateCarree(globe=cel_sphere)

    ra, dec = np.mgrid[-179.5:180:1, -89.5:90:1]
    fake_data_1 = generate_perlin_noise_2d(ra.shape, (1, 1))
    fake_data_2 = generate_perlin_noise_2d(ra.shape, (1, 1)) + 2

    norm = colors.Normalize(
        vmin=np.min([fake_data_1, fake_data_2]),
        vmax=np.max([fake_data_1, fake_data_2]),
    )

    for label, proj in [("ee", ccrs.EqualEarth),
                        ("mw", ccrs.Mollweide),
                        ("lc", ccrs.LambertCylindrical)]:
        for width in np.linspace(18, 22, 21):
            draw(fake_data_1, fake_data_2,
                 ra, dec, sky_plate, proj(globe=cel_sphere),
                 norm, width, 10, label)

def draw(d1, d2, ra, dec, data_crs, map_crs, norm, width, height, label):
    fig, (a1, a2) = plt.subplots(
        1, 2,
        figsize=(width, height),
        layout=ConstrainedLayoutEngine(
            h_pad=0, w_pad=0, hspace=0, wspace=0
        ),
        subplot_kw={
            "xlim": (-180, 180),
            "ylim": (-90, 90),
            "projection": map_crs,
        },
    )
    try:
        a1.gridlines(draw_labels=True)
        a2.gridlines(draw_labels=True)
        a1.contourf(ra, dec, d1,
                    transform=data_crs,
                    cmap="Greys",
                    norm=norm)
        a2.contourf(ra, dec, d2,
                    transform=data_crs,
                    cmap="Greys",
                    norm=norm)

        a1.set_title(label, loc="left")
        a2.set_title(f"{width}x{height}", loc="left")

        fig.colorbar(cm.ScalarMappable(norm=norm, cmap="Greys"),
                     shrink=0.5, pad=0.02, ax=[a1, a2])
        fig.savefig(f"layout_test_{label}_{width}x{height}.png",
                    bbox_inches="tight", pad_inches=0.125)
    finally:
        plt.close(fig)


# stolen from https://pvigier.github.io/2018/06/13/perlin-noise-numpy.html
def generate_perlin_noise_2d(shape, res):
    def f(t):
        return 6*t**5 - 15*t**4 + 10*t**3

    delta = (res[0] / shape[0], res[1] / shape[1])
    d = (shape[0] // res[0], shape[1] // res[1])
    grid = np.mgrid[0:res[0]:delta[0],0:res[1]:delta[1]].transpose(1, 2, 0) % 1
    # Gradients
    angles = 2*np.pi*np.random.rand(res[0]+1, res[1]+1)
    gradients = np.dstack((np.cos(angles), np.sin(angles)))
    g00 = gradients[0:-1,0:-1].repeat(d[0], 0).repeat(d[1], 1)
    g10 = gradients[1:,0:-1].repeat(d[0], 0).repeat(d[1], 1)
    g01 = gradients[0:-1,1:].repeat(d[0], 0).repeat(d[1], 1)
    g11 = gradients[1:,1:].repeat(d[0], 0).repeat(d[1], 1)
    # Ramps
    n00 = np.sum(grid * g00, 2)
    n10 = np.sum(np.dstack((grid[:,:,0]-1, grid[:,:,1])) * g10, 2)
    n01 = np.sum(np.dstack((grid[:,:,0], grid[:,:,1]-1)) * g01, 2)
    n11 = np.sum(np.dstack((grid[:,:,0]-1, grid[:,:,1]-1)) * g11, 2)
    # Interpolation
    t = f(grid)
    n0 = n00*(1-t[:,:,0]) + t[:,:,0]*n10
    n1 = n01*(1-t[:,:,0]) + t[:,:,0]*n11
    return np.sqrt(2)*((1-t[:,:,1])*n0 + t[:,:,1]*n1)


main()

Thanks to @SollyBunny for noticing that a large part of the problem was setting the xlim and ylim. These refer to the cartesian coordinates of the axes, which are distinct from the latitude and longitude. To visualise this, use plt.show() to see your figure in a gui window instead of saving to file. When you move your cursor over the gui window, you see two sets of numbers in the top right (if using QtAgg) or bottom right (if using TkAgg): these numbers are the x- and y- positions of the cursor, and the corresponding longitude and latitude. For EqualEarth, the map covers around ±155 in the x-direction and ±75 in the y-direction.

Just removing these xlim and ylim settings gets rid of the unwanted horizontal white space. You might also want to use ax.set_global() to ensure the full map is shown - without this the left and right of the map get slightly flattened around the equator.

            fig, ax = plt.subplots(
                figsize=(20, 10),
                layout=ConstrainedLayoutEngine(
                    h_pad=0, w_pad=0, hspace=0, wspace=0
                ),
                subplot_kw={"projection": proj(globe=cel_sphere)},
            )
            ax.set_global()

To remove the remaining unwanted vertical white space, you can save with bbox_inches='tight', which adjusts the figure size to fit around the artists on it.

fig.savefig(f"layout_test_{label}.png", bbox_inches="tight", pad_inches="layout")

With these changes I get