diff options
Diffstat (limited to 'export_mvt.py')
| -rw-r--r-- | export_mvt.py | 123 |
1 files changed, 111 insertions, 12 deletions
diff --git a/export_mvt.py b/export_mvt.py index d19909c..7b3137b 100644 --- a/export_mvt.py +++ b/export_mvt.py @@ -117,10 +117,10 @@ def createMVT(drv : gdal.Driver, path : str, return drv.Create(path, 0, 0, 0, **kwargs) # pylint: disable-next=too-many-branches -def exportSourceLayer(ds_src : gdal.Dataset, - lyr_src : ogr.Layer, +def exportSourceLayer(lyr_src : ogr.Layer, lyr_dst : ogr.Layer, layerdef : dict[str,Any], + fieldMap : tuple[list[str],list[int]], extent : ogr.Geometry|None = None) -> int: """Export a source layer.""" count0 = -1 @@ -168,7 +168,7 @@ def exportSourceLayer(ds_src : gdal.Dataset, spatialFilter = getSpatialFilterFromGeometry(extent, srs_src) transform_geometry = layerdef.get('transform-geometry', None) - columns = [ 'm.' + escape_identifier(lyr_src.GetFIDColumn()) ] + columns = [ 'm.' + escape_identifier(lyr_src.GetFIDColumn()) ] + fieldMap[0] geomFieldName_esc = escape_identifier(geomField.GetName()) if transform_geometry is None: columns.append('m.' + geomFieldName_esc) @@ -183,6 +183,7 @@ def exportSourceLayer(ds_src : gdal.Dataset, if cond is not None: query += ' WHERE ' + cond.strip() + ds_src = lyr_src.GetDataset() logging.debug('ExecuteSQL(%s%s)', query, '' if spatialFilter is None else ', spatialFilter=' + spatialFilter.ExportToWkt()) lyr_src = ds_src.ExecuteSQL(query, spatialFilter=spatialFilter) @@ -194,16 +195,26 @@ def exportSourceLayer(ds_src : gdal.Dataset, logging.debug('Source layer "%s" has %d features, of which %d are to be exported', layername, count0, count1) + fieldMap = fieldMap[1] + logging.debug('Field map: %s', str(fieldMap)) + + geom_type = lyr_src.GetGeomType() + bFlatten = geom_type == ogr.wkbUnknown or ogr.GT_HasM(geom_type) or ogr.GT_HasZ(geom_type) + bTransform = bFlatten or ct is not None + feature_count = 0 defn_dst = lyr_dst.GetLayerDefn() feature = lyr_src.GetNextFeature() while feature is not None: - geom = feature.GetGeometryRef().Clone() - if ct is not None and geom.Transform(ct) != ogr.OGRERR_NONE: - raise RuntimeError('Could not apply coordinate transformation') - geom.FlattenTo2D() feature2 = ogr.Feature(defn_dst) - feature2.SetGeometryDirectly(geom) + feature2.SetFromWithMap(feature, False, fieldMap) + if bTransform: + geom = feature2.GetGeometryRef() + if ct is not None and geom.Transform(ct) != ogr.OGRERR_NONE: + raise RuntimeError('Could not apply coordinate transformation') + if bFlatten: + geom.FlattenTo2D() + feature2.SetGeometryDirectly(geom) feature2.SetFID(feature.GetFID()) if lyr_dst.CreateFeature(feature2) != ogr.OGRERR_NONE: raise RuntimeError(f'Could not transfer source feature #{feature.GetFID()}') @@ -351,7 +362,84 @@ def exportMetadata(basedir : Path, data : dict[str,Any], finally: os.close(fd) -# pylint: disable-next=too-many-branches, too-many-statements +def getFieldMap(lyr_dst : ogr.Layer, lyr_src : ogr.Layer, + fieldMap : dict[str,str]|None) -> tuple[list[str],list[int]]: + """Create fields on the destination MVT layer, and return a list of + column statements along with a field map for the MVT export.""" + if fieldMap is None or len(fieldMap) == 0: + return [], [] + + if not lyr_dst.TestCapability(ogr.OLCCreateField): + raise RuntimeError(f'Destination layer "{lyr_dst.GetName()}" lacks ' + 'field creation capability') + + columns = {} + defn_src = lyr_src.GetLayerDefn() + drv_src = lyr_src.GetDataset().GetDriver() + for fld_dst, fld_src in fieldMap.items(): + idx_src = defn_src.GetFieldIndex(fld_src) + if idx_src < 0: + raise RuntimeError(f'Source layer "{lyr_src.GetName()}" has no field named "{fld_src}"') + + defn_dst = ogr.FieldDefn() + defn_src_fld = defn_src.GetFieldDefn(idx_src) + if fld_dst == 'ts': + if defn_src_fld.GetType() not in (ogr.OFTDate, ogr.OFTDateTime): + raise RuntimeError(f'Field "{fld_src}" of source layer "{lyr_src.GetName()}"' + ' has type ' + ogr.GetFieldTypeName(defn_src_fld.GetType()) + + ' (Date or DateTime expected)') + defn_dst.SetType(ogr.OFTInteger) + # signed int16 allows expressing dates from 1880-04-15 to 2059-09-18 + # which should be more than enough (it's not clear if the MVT format takes + # advantage of the reduced storage though) + defn_dst.SetSubType(ogr.OFSTInt16) + + if drv_src.ShortName == 'PostgreSQL': + column = 'CAST(m.' + escape_identifier(fld_src) + column += ' - date \'1970-01-01\' AS smallint)' + elif drv_src.ShortName in ('SQLite', 'GPKG'): + column = 'CAST(floor(julianday(m.' + escape_identifier(fld_src) + ')' + column += ' - 2440587.5) AS smallint)' + else: + raise NotImplementedError(f'Unsupported source driver {drv_src.ShortName} for ' + f'field "{fld_src}" (MVT field "{fld_dst}")') + + else: + raise NotImplementedError(f'Destination MVT field "{fld_dst}"') + + columns[fld_dst] = column + + defn_dst.SetName(fld_dst) + defn_dst.SetNullable(defn_src_fld.IsNullable()) + logging.debug('Create output field "%s" with type=%s, subtype=%s, nullable=%d', + defn_dst.GetName(), + ogr.GetFieldTypeName(defn_dst.GetType()), + ogr.GetFieldSubTypeName(defn_dst.GetSubType()), + defn_dst.IsNullable()) + + if lyr_dst.CreateField(defn_dst, approx_ok=False) != gdal.CE_None: + raise RuntimeError(f'Could not create field "{fld_dst}" ' + f'in destination MVT layer "{lyr_dst.GetName()}"') + + indices = {} + defn_dst = lyr_dst.GetLayerDefn() + for i in range(defn_dst.GetFieldCount()): + fld = defn_dst.GetFieldDefn(i) + name = fld.GetName() + if name in columns: + indices[name] = i + else: + logging.warning('Destination layer has unknown field #%d "%s"', i, name) + + ret = [None] * len(columns) + fieldMap = [-1] * defn_dst.GetFieldCount() + for idx, name in enumerate(columns.keys()): + i = indices[name] # intentionally crash if we didn't create that field + fieldMap[i] = idx + ret[idx] = columns[name] + ' AS ' + escape_identifier(name) + return (ret, fieldMap) + +# pylint: disable-next=too-many-branches, too-many-locals, too-many-statements def exportMVT(ds : gdal.Dataset, layers : dict[str,dict[str,Any]], sources : dict[str,Any], @@ -361,7 +449,8 @@ def exportMVT(ds : gdal.Dataset, drvname : str = 'MVT', default_options : dict[str,Any]|None = None, tile_extension : str = '.pbf', - compress : bool = False) -> None: + compress : bool = False, + compress_metadata : bool = False) -> None: """Export some layers to MVT.""" drv = gdal.GetDriverByName(drvname) if drv is None: @@ -369,6 +458,8 @@ def exportMVT(ds : gdal.Dataset, srs, extent = parseTilingScheme(default_options.get('tiling-scheme', None)) + last_modified_ns = max(last_modified.values()) * 1000000 if len(last_modified) > 0 else None + export_layers = {} mvtconf = {} for layername, layerdef in layers.items(): @@ -436,10 +527,14 @@ def exportMVT(ds : gdal.Dataset, if lyr_dst is None: raise RuntimeError(f'Could not create destination layer "{layername}"') + fieldMap = getFieldMap(lyr_dst, lyr_src, fieldMap=layerdef.get('fields', None)) + # TODO: GDAL exports features to a temporary SQLite database even though the source # is PostGIS hence is able to generate MVT with ST_AsMVT(). Letting PostGIS generate # tiles is likely to speed up things. - feature_count += exportSourceLayer(ds, lyr_src, lyr_dst, layerdef, extent=extent) + feature_count += exportSourceLayer(lyr_src, lyr_dst, layerdef, + fieldMap=fieldMap, + extent=extent) layer_count += 1 lyr_dst = None lyr_src = None @@ -494,7 +589,11 @@ def exportMVT(ds : gdal.Dataset, last_modified=last_modified, last_updated=creation_time // 1000000), dir_fd=dir_fd, - compress=compress) + compress=compress_metadata) + + if last_modified_ns is not None: + os.utime(mvtname, ns=(last_modified_ns, last_modified_ns), + dir_fd=dir_fd, follow_symlinks=False) try: # atomically exchange paths |