Skip to content

Instantly share code, notes, and snippets.

@attentive
Created October 13, 2020 21:11
Show Gist options
  • Save attentive/aaabcad2e25d729d35df4e4d4cc4f0ee to your computer and use it in GitHub Desktop.
Save attentive/aaabcad2e25d729d35df4e4d4cc4f0ee to your computer and use it in GitHub Desktop.
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
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 …")
@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
@attentive
Copy link
Author

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment