From 7cda119879cf48ba72ba34522fa9cdf9ef6d9b49 Mon Sep 17 00:00:00 2001 From: Guilhem Moulin Date: Sat, 21 Sep 2024 04:00:05 +0200 Subject: Add `webmap-publish` script to export layers to Mapbox Vector Tiles. --- common.py | 113 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 111 insertions(+), 2 deletions(-) (limited to 'common.py') 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 -- cgit v1.2.3