Add `webmap-publish` script to export layers to Mapbox Vector Tiles.

1 files changed, 111 insertions, 2 deletions

diff --git a/common.py b/common.py
index 79b5de8..1635195 100644
--- a/common.py
+++ b/common.py
@@ -65,7 +65,7 @@ def load_config(path=None, groupnames=None):
     for name, layerdefs in layers.items():
         if isinstance(layerdefs, dict) and 'sources' not in layerdefs:
             layers[name] = { 'sources': [layerdefs] }
-            for k in ['description', 'create']:
+            for k in ['description', 'create', 'publish']:
                 if k in layerdefs:
                     layers[name][k] = layerdefs.pop(k)
             layerdefs = layers[name]
@@ -190,7 +190,7 @@ def format_time(s):
 
 # Return a boolean indicating whether the installer GDAL version is
 # greater than or equal to the provider (maj, min, rev) triplet.
-def gdal_version_min(maj=0, min=0, rev=0):
+def gdalVersionMin(maj=0, min=0, rev=0):
     if maj < 1 or (maj == 1 and min < 10):
         # GDAL_VERSION_NUM() macro was changed in 1.10.  That version
         # was released in 2013 so we blindly assume the installer
@@ -198,11 +198,120 @@ def gdal_version_min(maj=0, min=0, rev=0):
         return True
 
     from osgeo import gdal
+    gdal.UseExceptions()
+
     version_cur = int(gdal.VersionInfo());
     # cf. GDAL_COMPUTE_VERSION(maj,min,rev) in gcore/gdal_version.h.in
     version_min = maj*1000000 + min*10000 + rev*100
     return version_min <= version_cur
 
+# Wrapper around gdal.MajorObject.GetMetadataItem(name)
+def gdalGetMetadataItem(o, k):
+    v = o.GetMetadataItem(k)
+    if v is not None and isinstance(v, str):
+        return v.upper() == 'YES'
+    else:
+        return False
+
+# Escape the given identifier, cf.
+# swig/python/gdal-utils/osgeo_utils/samples/validate_gpkg.py:_esc_id()
+def escapeIdentifier(identifier):
+    if '\x00' in identifier:
+        raise Exception(f'Invalid identifier "{identifier}"')
+    # SQL:1999 delimited identifier
+    return '"' + identifier.replace('"', '""') + '"'
+
+# Return a pair kwargs and driver to use with gdal.OpenEx()
+def gdalSetOpenExArgs(option_dict, flags=0):
+    from osgeo import gdal
+    gdal.UseExceptions()
+
+    kwargs = { 'nOpenFlags': gdal.OF_VECTOR | flags }
+
+    fmt = option_dict.get('format', None)
+    if fmt is None:
+        drv = None
+    else:
+        drv = gdal.GetDriverByName(fmt)
+        if drv is None:
+            raise Exception(f'Unknown driver name "{fmt}"')
+        elif not gdalGetMetadataItem(drv, gdal.DCAP_VECTOR):
+            raise Exception(f'Driver "{drv.ShortName}" has no vector capabilities')
+        kwargs['allowed_drivers'] = [ drv.ShortName ]
+
+    oo = option_dict.get('open-options', None)
+    if oo is not None:
+        kwargs['open_options'] = [ k + '=' + str(v) for k, v in oo.items() ]
+    return kwargs, drv
+
+# Return the decoded Spatial Reference System
+def getSRS(srs_str):
+    if srs_str is None:
+        return
+
+    from osgeo import osr
+    osr.UseExceptions()
+
+    srs = osr.SpatialReference()
+    if srs_str.startswith('EPSG:'):
+        code = int(srs_str.removeprefix('EPSG:'))
+        srs.ImportFromEPSG(code)
+    else:
+        raise Exception(f'Unknown SRS {srs_str}')
+    logging.debug('Default SRS: "%s" (%s)', srs.ExportToProj4(), srs.GetName())
+    return srs
+
+# Convert extent [minX, minY, maxX, maxY] into a polygon and assign the
+# given SRS.  Return a pair with the densified and non-densified extent.
+# Like apps/ogr2ogr_lib.cpp, the former is obtained by segmentizing the
+# polygon to make sure it is sufficiently densified when transforming to
+# source layer SRS for spatial filtering.
+def getExtent(extent, srs=None):
+    if extent is None:
+        return None, None
+
+    if not (isinstance(extent, list) or isinstance(extent, tuple)) or len(extent) != 4:
+        raise Exception(f'Invalid extent {extent}')
+    elif srs is None:
+        raise Exception('Configured extent but no SRS')
+
+    logging.debug('Configured extent in %s: %s',
+        srs.GetName(), ', '.join(map(str, extent)))
+
+    from osgeo import ogr, osr
+    ogr.UseExceptions()
+
+    ring = ogr.Geometry(ogr.wkbLinearRing)
+    ring.AddPoint_2D(extent[0], extent[1])
+    ring.AddPoint_2D(extent[2], extent[1])
+    ring.AddPoint_2D(extent[2], extent[3])
+    ring.AddPoint_2D(extent[0], extent[3])
+    ring.AddPoint_2D(extent[0], extent[1])
+
+    polygon = ogr.Geometry(ogr.wkbPolygon)
+    polygon.AddGeometry(ring)
+
+    # we expressed extent as minX, minY, maxX, maxY (easting/northing
+    # ordered, i.e., in traditional GIS order)
+    srs2 = srs.Clone()
+    srs2.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    polygon.AssignSpatialReference(srs2)
+    if not srs2.IsSame(srs):
+        polygon.TransformTo(srs)
+
+    # densify the rectangle to avoid issues when reprojecting to the
+    # source layer SRS, cf. apps/ogr2ogr_lib.cpp:ApplySpatialFilter()
+    polygon_dense = polygon.Clone()
+    segment_distance_metre = 10 * 1000
+    if srs.IsGeographic():
+        dfMaxLength = segment_distance_metre / math.radians(srs.GetSemiMajor())
+        polygon_dense.Segmentize(dfMaxLength)
+    elif srs.IsProjected():
+        dfMaxLength = segment_distance_metre / srs.GetLinearUnits()
+        polygon_dense.Segmentize(dfMaxLength)
+
+    return polygon_dense, polygon
+
 
 ######
 # The function definitions below are taken from cpython's source code