Using GDAL and Python (QGIS 3.10 compatible) to carve a set of MBTiles from a larger raster based on the polygon region features in a shapefile

carve_mbtiles.py

import itertools, os, time
from sys import stdout
from osgeo import gdal, ogr, osr

# change these to the local paths if running this yourself
scratch = r"D:\scratch"
mbtilesFile = r"D:\big-mbtiles-file.mbtiles"
shapefile = r"C:\temp\regions.shp"

webMercator = osr.SpatialReference()
webMercator.ImportFromEPSG(3857)

driver = ogr.GetDriverByName("ESRI Shapefile")
regions = driver.Open(shapefile, 0) 

# check to see if shapefile is found
if regions is None:
    print(f"Could not open regions {shapefile}")
else:
    print(f"Opened {shapefile}")

# check to see if MBtiles is found
try: 
    mbtiles = gdal.Open(mbtilesFile)
except RuntimeError as e:
    print(f"Could not open MBTiles input file {mbtilesFile}")
    print(str(e))

layer = regions.GetLayer()
layerSrs = layer.GetSpatialRef()
transform = osr.CoordinateTransformation(layerSrs, webMercator)
regionCount = layer.GetFeatureCount()

start = time.process_time()

def progressCallback(complete, message, data):
    print(f"{(complete * 100):.2f}%, {(time.process_time() - start):.2f}s", end='\r')
    stdout.flush()
    return 1

for i, region in enumerate(layer):
    regionId = region.GetField("AREA_ID")
    print(f"Carving region '{regionId}' ({i + 1} of {regionCount}) …")

    # get region bounds in Web Mercator
    geom = region.GetGeometryRef()
    geom.Transform(transform)
    envelope = geom.GetEnvelope()
    bounds = f"{envelope[0]}, {envelope[2]}, {envelope[1]}, {envelope[3]}"
    print(f"Bounds: {bounds}")
    
    vrtOptions = {
        "outputBounds": (envelope[0], envelope[2], envelope[1], envelope[3]),
        "outputSRS": "EPSG:3857"
    }
    
    # create region-windowed input VRT    
    vrtFile = os.path.join(scratch, f"{regionId}.vrt")
    vrt = gdal.BuildVRT(vrtFile, [mbtilesFile], options=gdal.BuildVRTOptions(**vrtOptions))
    vrt.FlushCache()
    
    # create initial translated MBTiles file
    regionMbtilesFile = os.path.join(scratch, f"{regionId}.mbtiles")
    if not os.path.exists(regionMbtilesFile):
        print(f"Creating MBTiles for '{regionId}' …")
        print(f"Output file: {regionMbtilesFile}")
        start = time.process_time()
        translateOptions = {
            "callback": progressCallback,
            "format": "MBTiles",
            "creationOptions": ["MINZOOM=0","TILE_FORMAT=JPEG"]
        }
        gdal.Translate(regionMbtilesFile, vrt,
                       options=gdal.TranslateOptions(**translateOptions))
    
        print("\nMBTiles created.")
    
        # add overviews
        print("Adding overviews …")
        regionMbtiles = gdal.Open(regionMbtilesFile, gdal.GA_Update)
        overviewList = [2**i for i in range(1, 13)]
        regionMbtiles.BuildOverviews("AVERAGE", overviewList, callback=progressCallback)
        print("\nOverviews added.")

        # close dataset
        regionMbtiles = None
    else:
        print(f"Skipping MBTiles creation for '{regionId}', already exists …")

setup_carve_mbtiles.bat

@ECHO OFF

::Change this to your QGIS 3.x version …
set OSGEO4W_ROOT=C:\Program Files\QGIS 3.10

set PATH=%OSGEO4W_ROOT%\bin;%PATH%
set PATH=%PATH%;%OSGEO4W_ROOT%\apps\qgis\bin

@echo off
call "%OSGEO4W_ROOT%\bin\o4w_env.bat"
call "%OSGEO4W_ROOT%\bin\qt5_env.bat"
call "%OSGEO4W_ROOT%\bin\py3_env.bat"
@echo off
path %OSGEO4W_ROOT%\apps\qgis-ltr\bin;%OSGEO4W_ROOT%\apps\grass\grass78\lib;%OSGEO4W_ROOT%\apps\grass\grass78\bin;%PATH%

cd /d %~dp0

This should be compatible with QGIS 3.x on Windows 10, just install it, fire up a command prompt and run the batch file to set up your environment, then you can use python carve_mbtiles.py to run the Python if you like.

Same code will work fine on other platforms with minor modifications.