Lambert projection with region of interest

I’m trying to make a equal-area map over a few ROIs with minimal local distortion. It looks like Lambert azimuthal should work, but the result always shows the entire globe, not just the ROI. Here’s an example taken from the tutorial on regions, just changing the projection to A15/41/15c. I’d like the map to show only what is inside the box. Do you have any suggestion?

fig.coast(
    region="10/20/35/45",
    projection="A15/41/15c",
    land="lightgray",
    water="white",
    borders="1/0.5p",
    shorelines="1/0.5p",
    frame="ag",
)
fig.plot(x=[10,20,20,10,10], y=[35,35,45,45,35], pen="2p,black")
fig.show()
fig = pygmt.Figure()

You can try these two approaches:

(I)
You can set the horizon via the projection parameter to limit the displayed area, please see
Azimuthal equidistant projection — PyGMT.

import pygmt

fig = pygmt.Figure()

# Set the horizon to 10 deg
fig.basemap(region=[10, 20, 35, 45], projection="A15/40/10/15c", frame="afg")
fig.coast(shorelines="1/0.5p,black")

fig.plot(x=[10, 20, 20, 10, 10], y=[35, 35, 45, 45, 35], pen="3p,blue", no_clip=True)

fig.show()

(II)
You can make your map have a rectangular shape using +r with the region parameter, please see 6. GMT Map Projections — GMT 6.6.0 documentation

import pygmt

fig = pygmt.Figure()

# Use +r to get a rectangular map
# Note region has to be given as the coordinates of the lower left and the
#      upper right corners
fig.basemap(region="10/35/20/45+r", projection="A15/40/15c", frame="afg")
fig.coast(shorelines="1/0.5p,black")

fig.plot(x=[10, 20, 20, 10, 10], y=[35, 35, 45, 45, 35], pen="3p,blue", no_clip=True)

fig.show()