Putting subplot-titles in Density Mapbox charts

So, I am trying to make a chart using subplots of go.Densitymapbox however I am having the issue that I want to title the subplots with month names. My code looks like:

gear = 'trawlers'
flag = "GBR"

df_UK_fishing = df[(df.flag == flag) & (df.geartype == gear)]
df_UK_fishing['flag'].value_counts()

data = []
months = df_UK_fishing['Month'].unique()
# months = ['January', 'February']

fig = go.Figure()

COLS = 4
ROWS = 3
x = 0
y = ROWS
for i in range(len(months)):
    df_tmp = df_UK_fishing[df_UK_fishing['Month'] == months[i]]
    mapbox = 'mapbox'+str(i+1) if i != 0 else 'mapbox'
    fig.add_trace(
        go.Densitymapbox(
            lat=df_tmp.cell_ll_lat, 
            lon=df_tmp.cell_ll_lon,
            name = months[i],
            z=df_tmp.fishing_hours, 
            radius = 4,
            subplot=mapbox,
            zmin=0,
            zmax=25,
        )
    )
    fig.update_layout({
        mapbox:{ 
            'style': 'open-street-map',
            'domain': {'x': [float(x)/float(COLS), float(x+1)/float(COLS)], 'y': [float(y-1)/float(ROWS), float(y)/float(ROWS)]},
            'center': {'lat': cen_lat, 'lon':cen_lon},
            'zoom': 4,
#             'layers': [{
#                 "sourcetype": "raster",
#                 "sourceattribution": str(months[i]),
#                 "source": [
#                     "https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}"
#                 ]
#             }]
        },
        'meta':{
            mapbox:{
                'title':{'text':months[i]},
                'dragmode': False,
                'margin': {"r":30, "t":30,"l":30,"b":30},
                'showlegend': True,
            }
        }
    })
    
    x = x+1
    if x >= COLS:
        y=y-1
        x = 0
    if y <= 0:
        y = ROWS

fig.update_layout(
    margin={"r":0, "t":30,"l":0,"b":0},
    title = f"2019 monthly fishing data for {flag} {gear}", 
)

Which creates a chart that looks like this:

It looks good and I donโ€™t mind the โ€œOpen Street Mapsโ€ mark but what Iโ€™d really like, recreated by changing 'style': 'white-bg' and uncommenting the 'layers' from fig.update_layout(), is to produce a chart like this:

Basically what I want is to have Months as the sub-plot titles.

Any hot tips on to get that working?

I figured it out that I hadnโ€™t set the layers traces to โ€˜belowโ€™ so the data is there but the map is printed above the data.

The code looks like:

def chart_fishing(gear, flag, is_flag = True):
    if is_flag:
        df_UK_fishing = df[(df.flag == flag) & (df.geartype == gear)]
    else:
        df_UK_fishing = df[(df.flag != flag) & (df.geartype == gear)]
        flag = "not " + flag
    
    gear = gear.replace("_", " ").title()
    
    if not df_UK_fishing.empty:
        data = []
        months = df_UK_fishing['Month'].unique()

        fig = go.Figure()

        COLS = 4
        ROWS = 3
        x = 0
        y = ROWS
        for i in range(len(months)):
            df_tmp = df_UK_fishing[df_UK_fishing['Month'] == months[i]]
            mapbox = 'mapbox'+str(i+1) if i != 0 else 'mapbox'

            fig.add_trace(
                go.Densitymapbox(
                    lat=df_tmp.cell_ll_lat, 
                    lon=df_tmp.cell_ll_lon,
                    z=df_tmp.fishing_hours, 
                    radius = 4,
                    subplot=mapbox,
                    zmin=0,
                    zmax=25,
#                     colorscale='viridis',
                )
            )

            fig.update_layout({
                mapbox:{ 
                    'domain': {
                        'x': [float(x)/float(COLS), float(x+1)/float(COLS)], 
                        'y': [float(y-1)/float(ROWS), float(y)/float(ROWS)]
                    },
                    'center': {'lat': cen_lat, 'lon':cen_lon},
                    'zoom': 4,
                    'style': 'white-bg',
                    'layers': [{
                        'below': 'traces', # This is the bit of code I forgot
                        "sourcetype": "raster",
                        "sourceattribution": str(months[i]),
                        "source": [
"https://basemap.nationalmap.gov/arcgis/rest/services/USGSHydroCached/MapServer/tile/{z}/{y}/{x}"
                        ]
                    }]
                },
            })

            x = x+1
            if x >= COLS:
                y=y-1
                x = 0
            if y <= 0:
                y = ROWS
        
        nations = df_UK_fishing.Nationality.unique().tolist()
        title_height = 0
        for i, nat in enumerate(nations):
            if len(nations) > 3 and ((i+1)%8 == 0 or len(nations[i]) > 50):
                nations[i] = "<br>" + nations[i]
                title_height = title_height + 1 
                
        fig.update_layout(
            margin={"r":0, "t":60+30*title_height,"l":0,"b":0},
            title = f"2019 monthly fishing hours for {gear} <br>{', '.join(nations)}",
            title_font_size = 15,
        )
        im_file = f"D:\Data\Fishing Data Science\Data\Maps\\2019_{flag}_{gear}.png"
        pio.write_image(fig, im_file, format='png', scale=5)

        img = mpimg.imread(im_file)
        plt.figure(figsize=(20,20))
        plt.imshow(img)
#         fig.show()
    else:
        print("No data")

chart_fishing('trawlers', 'GBR', False)

Now my charts look like this: