Add `webmap-import` script to import source layers.

There is still a few things to do (such as reprojection and geometry
changes) but it's mostly working.

We roll out our own ogr2ogr/GDALVectorTranslate()-like function version
because GDALVectorTranslate() insists in calling StartTransaction()
https://github.com/OSGeo/gdal/issues/3403 while we want a single
transaction for the entire desination layer, including truncation,
source imports, and metadata changes.

Surprisingly our version is not much slower than the C++ one.  Importing
the 157446 (of 667034) features from sksUtfordAvverk-2000-2015.shp takes
14.3s while

	ogr2ogr -f PostgreSQL \
	    -doo ACTIVE_SCHEMA=postgis \
	    --config PG_USE_COPY YES \
	    --config OGR_TRUNCATE YES \
	    -append \
	    -fieldmap "0,-1,-1,-1,-1,1,2,3,4,5,6,7,8,9,10,11,12,13" \
	    -nlt MULTIPOLYGON -nlt PROMOTE_TO_MULTI \
	    -gt unlimited \
	    -spat 110720 6927136 1159296 7975712 \
	    -nln "sks:UtfordAvverk" \
	    PG:"dbname='webmap' user='webmap_import'" \
	    /tmp/x/sksUtfordAvverk-2000-2015.shp \
	    sksUtfordAvverk-2000-2015

takes 14s.

Merely opening /tmp/x/sksUtfordAvverk-2000-2015.shp and looping through
its (extent-filtered) features results in a runtime of 4.3s.

1 files changed, 1017 insertions, 0 deletions

diff --git a/webmap-import b/webmap-import
new file mode 100755
index 0000000..6daceee
--- /dev/null
+++ b/webmap-import
@@ -0,0 +1,1017 @@
+#!/usr/bin/python3
+
+#----------------------------------------------------------------------
+# Backend utilities for the Klimatanalys Norr project (extract/import layers)
+# Copyright © 2024 Guilhem Moulin <info@guilhem.se>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+#----------------------------------------------------------------------
+
+import os
+import logging
+import argparse
+import tempfile
+from fnmatch import fnmatchcase
+from pathlib import Path
+
+from osgeo import gdal, ogr, osr
+from osgeo.gdalconst import (
+    OF_VECTOR as GDAL_OF_VECTOR,
+    OF_ALL as GDAL_OF_ALL,
+    OF_READONLY as GDAL_OF_READONLY,
+    OF_UPDATE as GDAL_OF_UPDATE,
+    OF_VERBOSE_ERROR as GDAL_OF_VERBOSE_ERROR,
+    CE_None as GDAL_CE_None,
+    DCAP_CREATE as GDAL_DCAP_CREATE,
+    DCAP_VECTOR as GDAL_DCAP_VECTOR,
+    DCAP_DEFAULT_FIELDS as GDAL_DCAP_DEFAULT_FIELDS,
+    DCAP_NOTNULL_FIELDS as GDAL_DCAP_NOTNULL_FIELDS,
+    DCAP_UNIQUE_FIELDS as GDAL_DCAP_UNIQUE_FIELDS,
+)
+import osgeo.gdalconst as gdalconst
+gdal.UseExceptions()
+
+import common
+
+# Wrapper around gdal.MajorObject.GetMetadataItem(name)
+def getMetadataItem(o, k):
+    v = o.GetMetadataItem(k)
+    if v is not None and isinstance(v, str):
+        return v.upper() == 'YES'
+    else:
+        return False
+
+# Return kwargs and driver for OpenEx()
+def setOpenExArgs(option_dict, flags=0):
+    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 getMetadataItem(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
+
+# Open and return the output DS.  It is created if create=False or
+# create-options is a non-empty dictionary.
+def openOutputDS(def_dict):
+    path = def_dict['path']
+    kwargs, drv = setOpenExArgs(def_dict, flags=GDAL_OF_UPDATE|GDAL_OF_VERBOSE_ERROR)
+    try:
+        logging.debug('OpenEx(%s, %s)', path, str(kwargs))
+        return gdal.OpenEx(path, **kwargs)
+    except RuntimeError as e:
+        if not (gdal.GetLastErrorType() >= gdalconst.CE_Failure and
+                gdal.GetLastErrorNo() == gdalconst.CPLE_OpenFailed):
+            # not an open error
+            raise e
+
+        dso2 = None
+        try:
+            dso2 = gdal.OpenEx(path, nOpenFlags=GDAL_OF_ALL | GDAL_OF_UPDATE)
+        except Exception:
+            pass
+        if dso2 is not None:
+            # path exists but can't be open with OpenEx(path, **kwargs)
+            raise e
+
+        try:
+            dso2 = gdal.OpenEx(path, nOpenFlags=GDAL_OF_ALL)
+        except Exception:
+            pass
+        if dso2 is not None:
+            # path exists but can't be open with OpenEx(path, **kwargs)
+            raise e
+
+        dsco = def_dict.get('create-options', None)
+        if not def_dict.get('create', False) and dsco is None:
+            # not configured for creation
+            raise e
+        if drv is None or not getMetadataItem(drv, GDAL_DCAP_CREATE):
+            # not capable of creation
+            raise e
+
+    if 'open_options' in kwargs:
+        # like ogr2ogr(1)
+        logging.warning('Destination\'s open options ignored ' +
+            'when creating the output datasource')
+
+    kwargs2 = { 'eType': gdal.GDT_Unknown }
+    if dsco is not None:
+        kwargs2['options'] = [ k + '=' + str(v) for k, v in dsco.items() ]
+
+    logging.debug('Create(%s, %s, eType=%s%s)', drv.ShortName, path, kwargs2['eType'],
+        ', options=' + str(kwargs2['options']) if 'options' in kwargs2 else '')
+    # XXX racy, a GDAL equivalent of O_EXCL would be nice
+    return drv.Create(path, 0, 0, 0, **kwargs2)
+
+# cf. ogr/ogrgeometry.cpp:OGRFromOGCGeomType()
+def fromGeomTypeName(name):
+    name = name.upper()
+
+    isMeasured = False
+    if name.endswith('M'):
+        isMeasured = True
+        name = name.removesuffix('M')
+
+    convertTo3D = False
+    if name.endswith('Z'):
+        convertTo3D = True
+        name = name.removesuffix('Z')
+
+    if name == 'POINT':
+        eGType = ogr.wkbPoint
+    elif name == 'LINESTRING':
+        eGType = ogr.wkbLineString
+    elif name == 'POLYGON':
+        eGType = ogr.wkbPolygon
+    elif name == 'MULTIPOINT':
+        eGType = ogr.wkbMultiPoint
+    elif name == 'MULTILINESTRING':
+        eGType = ogr.wkbMultiLineString
+    elif name == 'MULTIPOLYGON':
+        eGType = ogr.wkbMultiPolygon
+    elif name == 'GEOMETRYCOLLECTION':
+        eGType = ogr.wkbGeometryCollection
+    elif name == 'CIRCULARSTRING':
+        eGType = ogr.wkbCircularString
+    elif name == 'COMPOUNDCURVE':
+        eGType = ogr.wkbCompoundCurve
+    elif name == 'CURVEPOLYGON':
+        eGType = ogr.wkbCurvePolygon
+    elif name == 'MULTICURVE':
+        eGType = ogr.wkbMultiCurve
+    elif name == 'MULTISURFACE':
+        eGType = ogr.wkbMultiSurface
+    elif name == 'TRIANGLE':
+        eGType = ogr.wkbTriangle
+    elif name == 'POLYHEDRALSURFACE':
+        eGType = ogr.wkbPolyhedralSurface
+    elif name == 'TIN':
+        eGType = ogr.wkbTIN
+    elif name == 'CURVE':
+        eGType = ogr.wkbCurve
+    elif name == 'SURFACE':
+        eGType = ogr.wkbSurface
+    else:
+        eGType = ogr.wkbUnknown
+
+    if convertTo3D:
+        eGType = ogr.GT_SetZ(eGType)
+
+    if isMeasured:
+        eGType = ogr.GT_SetM(eGType)
+
+    return eGType
+
+# Parse geometry type, cf. ogr2ogr_lib.cpp
+def parseGeomType(name):
+    if name is None:
+        return ogr.wkbUnknown
+    name2 = name.upper()
+
+    is3D = False
+    if name2.endswith('25D'):
+        name2 = name2[:-3] # alias
+        is3D = True
+    elif name2.endswith('Z'):
+        name2 = name2[:-1]
+        is3D = True
+
+    if name2 == 'NONE':
+        eGType = ogr.wkbNone
+    elif name2 == 'GEOMETRY':
+        eGType = ogr.wkbUnknown
+    else:
+        eGType = fromGeomTypeName(name2)
+        if eGType == ogr.wkbUnknown:
+            raise Exception(f'Unknown geometry type "{name}"')
+
+    if eGType != ogr.wkbNone and is3D:
+        eGType = ogr.GT_SetZ(eGType)
+
+    return eGType
+
+# cf. ogr/ogr_core.h's enum OGRFieldType;
+def getFieldTypeCode(name):
+    if name is None:
+        raise Exception('getFieldTypeCode(None)')
+    name2 = name.lower()
+    if name2 == 'integer':
+        # simple 32bit integer
+        return ogr.OFTInteger
+    elif name2 == 'integerlist':
+        # List of 32bit integers
+        return ogr.OFTIntegerList
+    elif name2 == 'real':
+        # Double Precision floating point
+        return ogr.OFTReal
+    elif name2 == 'reallist':
+        # List of doubles
+        return ogr.OFTRealList
+    elif name2 == 'string':
+        # String of ASCII chars
+        return ogr.OFTString
+    elif name2 == 'stringlist':
+        # Array of strings
+        return ogr.OFTStringList
+    elif name2 == 'binary':
+        # Raw Binary data
+        return ogr.OFTBinary
+    elif name2 == 'date':
+        # Date
+        return ogr.OFTDate
+    elif name2 == 'time':
+        # Time
+        return ogr.OFTTime
+    elif name2 == 'datetime':
+        # Date and Time
+        return ogr.OFTDateTime
+    elif name2 == 'integer64':
+        # Single 64bit integer
+        return ogr.OFTInteger64
+    elif name2 == 'integer64list':
+        # List of 64bit integers
+        return ogr.OFTInteger64List
+    else:
+        raise Exception(f'Unknown field type "{name}"')
+
+# cf. ogr/ogr_core.h's enum OGRFieldSubType;
+def getFieldSubTypeCode(name):
+    if name is None:
+        raise Exception('getFieldSubTypeCode(None)')
+    name2 = name.lower()
+    if name2 == 'none':
+        # No subtype. This is the default value.
+        return ogr.OFSTNone
+    elif name2 == 'bool':
+        # Boolean integer. Only valid for OFTInteger and OFTIntegerList.
+        return ogr.OFSTBoolean
+    elif name2 == 'int16':
+        # Signed 16-bit integer. Only valid for OFTInteger and OFTIntegerList.
+        return ogr.OFSTInt16
+    elif name2 == 'float32':
+        # Single precision (32 bit) floating point. Only valid for OFTReal and OFTRealList.
+        return ogr.OFSTFloat32
+    elif name2 == 'json':
+        # JSON content. Only valid for OFTString.
+        return ogr.OFSTJSON
+    elif name2 == 'uuid':
+        # UUID string representation. Only valid for OFTString.
+        return ogr.OFSTUUID
+    else:
+        raise Exception(f'Unknown field subtype "{name2}"')
+
+# Validate layer creation options and schema.  The schema is modified in
+# place with the parsed result.
+# (We need the driver of the output dataset to determine capability on
+# constraints.)
+def validateSchema(layers, drvo=None, lco_defaults=None):
+    # cache driver capabilities
+    drvoSupportsDefaultFields = getMetadataItem(drvo, GDAL_DCAP_DEFAULT_FIELDS)
+    drvoSupportsNotNULLFields = getMetadataItem(drvo, GDAL_DCAP_NOTNULL_FIELDS)
+    drvoSupportsUniqueFields  = getMetadataItem(drvo, GDAL_DCAP_UNIQUE_FIELDS)
+
+    for layername, layerdef in layers.items():
+        create = layerdef.get('create', None)
+        if create is None or len(create) < 1:
+            logging.warning('Layer "%s" has no creation schema', layername)
+            continue
+
+        # prepend global layer creation options (dataset:create-layer-options)
+        # and build the option=value list
+        lco = create.get('options', None)
+        if lco_defaults is not None or lco is not None:
+            options = []
+            if lco_defaults is not None:
+                options += [ k + '=' + str(v) for k, v in lco_defaults.items() ]
+            if lco is not None:
+                options += [ k + '=' + str(v) for k, v in lco.items() ]
+            create['options'] = options
+
+        # parse geometry type
+        create['geometry-type'] = parseGeomType(create.get('geometry-type', None))
+
+        fields = create.get('fields', None)
+        if fields is None:
+            create['fields'] = []
+        else:
+            fields_set = set()
+            for idx, fld_def in enumerate(fields):
+                fld_name = fld_def.get('name', None)
+                if fld_name is None or fld_name == '':
+                    raise Exception(f'Field #{idx} has no name')
+                if fld_name in fields_set:
+                    raise Exception(f'Duplicate field "{fld_name}"')
+                fields_set.add(fld_name)
+
+                fld_def2 = { 'Name': fld_name }
+                for k, v in fld_def.items():
+                    k2 = k.lower()
+                    if k2 == 'name':
+                        pass
+                    elif k2 == 'alternativename' or k2 == 'alias':
+                        fld_def2['AlternativeName'] = v
+                    elif k2 == 'comment':
+                        fld_def2['Comment'] = v
+
+                    elif k2 == 'type':
+                        fld_def2['Type'] = getFieldTypeCode(v)
+                    elif k2 == 'subtype':
+                        fld_def2['SubType'] = getFieldSubTypeCode(v)
+                    elif k2 == 'tzflag':
+                        pass # TODO
+                        #fld_def2['TZFlag'] = v
+                    elif k2 == 'width' and v is not None and isinstance(v, int):
+                        fld_def2['Width'] = v
+                    elif k2 == 'precision' and v is not None and isinstance(v, int):
+                        fld_def2['Precision'] = v
+
+                    # constraints
+                    elif k2 == 'default':
+                        if drvoSupportsDefaultFields:
+                            fld_def2['Default'] = v
+                        else:
+                            logging.warning('%s driver lacks GDAL_DCAP_DEFAULT_FIELDS support',
+                                drvo.ShortName)
+                    elif k2 == 'nullable' and v is not None and isinstance(v, bool):
+                        if drvoSupportsNotNULLFields:
+                            fld_def2['Nullable'] = v
+                        else:
+                            logging.warning('%s driver lacks GDAL_DCAP_NOTNULL_FIELDS support',
+                                drvo.ShortName)
+                    elif k2 == 'unique' and v is not None and isinstance(v, bool):
+                        if drvoSupportsUniqueFields:
+                            fld_def2['Unique'] = v
+                        else:
+                            logging.warning('%s driver lacks GDAL_DCAP_UNIQUE_FIELDS support',
+                                drvo.ShortName)
+                    else:
+                        raise Exception(f'Field "{fld_name}" has unknown key "{k}"')
+
+                fields[idx] = fld_def2
+
+# Return the decoded Spatial Reference System
+def getSRS(srs_str):
+    if srs_str is None:
+        return
+    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
+
+# Validate the output layer against the provided SRS and creation options
+def validateOutputLayer(lyr, srs=None, options=None):
+    ok = True
+
+    # ensure the output SRS is equivalent
+    if srs is not None:
+        srs2 = lyr.GetSpatialRef()
+        # cf. apps/ogr2ogr_lib.cpp
+        srs_options = [
+            'IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES',
+            'CRITERION=EQUIVALENT'
+        ]
+        if not srs.IsSame(srs2, srs_options):
+            logging.warning('Output layer "%s" has SRS %s,\nexpected %s',
+                lyr.GetName(),
+                srs2.ExportToPrettyWkt(),
+                srs.ExportToPrettyWkt())
+            ok = False
+
+    if options is None:
+        return ok
+
+    layerDefn = lyr.GetLayerDefn()
+    n = layerDefn.GetGeomFieldCount()
+    if n != 1:
+        logging.warning('Output layer "%s" has %d != 1 geometry fields', layername, n)
+
+    geom_type1 = lyr.GetGeomType()
+    geom_type2 = options['geometry-type']
+    if geom_type1 != geom_type2:
+        logging.warning('Output layer "%s" has geometry type #%d (%s), expected #%d (%s)',
+            lyr.GetName(),
+            geom_type1, ogr.GeometryTypeToName(geom_type1),
+            geom_type2, ogr.GeometryTypeToName(geom_type2))
+        ok = False
+
+    fields = options.get('fields', None)
+    if fields is not None:
+        for fld in fields:
+            fldName = fld['Name']
+
+            idx = layerDefn.GetFieldIndex(fldName)
+            if idx < 0:
+                logging.warning('Output layer "%s" has no field named "%s"',
+                    lyr.GetName(), fldName)
+                ok = False
+                continue
+            defn = layerDefn.GetFieldDefn(idx)
+
+            if 'AlternativeName' in fld:
+                v1 = defn.GetAlternativeName()
+                v2 = fld['AlternativeName']
+                if v1 != v2:
+                    logging.warning('Field "%s" has AlternativeName="%s", expected "%s"',
+                        fldName, v1, v2)
+                    ok = False
+
+            if 'Comment' in fld:
+                v1 = defn.GetComment()
+                v2 = fld['Comment']
+                if v1 != v2:
+                    logging.warning('Field "%s" has Comment="%s", expected "%s"',
+                        fldName, v1, v2)
+                    ok = False
+
+            if 'Type' in fld:
+                v1 = defn.GetType()
+                v2 = fld['Type']
+                if v1 != v2:
+                    logging.warning('Field "%s" has Type=%d (%s), expected %d (%s)',
+                        fldName,
+                        v1, ogr.GetFieldTypeName(v1),
+                        v2, ogr.GetFieldTypeName(v2))
+                    ok = False
+
+            if 'SubType' in fld:
+                v1 = defn.GetSubType()
+                v2 = fld['SubType']
+                if v1 != v2:
+                    logging.warning('Field "%s" has SubType=%d (%s), expected %d (%s)',
+                        fldName,
+                        v1, ogr.GetFieldSubTypeName(v1),
+                        v2, ogr.GetFieldSubTypeName(v2))
+                    ok = False
+
+            if 'TZFlag' in fld:
+                v1 = defn.GetTZFlag()
+                v2, n2 = fld['TZFlag']
+                if v1 != v2:
+                    logging.warning('Field "%s" has TZFlag=%d, expected %d (%s)',
+                        fldName, v1, v2, n2)
+                    ok = False
+
+            if 'Precision' in fld:
+                v1 = defn.GetPrecision()
+                v2 = fld['Precision']
+                if v1 != v2:
+                    logging.warning('Field "%s" has Precision=%d, expected %d',
+                        fldName, v1, v2)
+                    ok = False
+
+            if 'Width' in fld:
+                v1 = defn.GetWidth()
+                v2 = fld['Width']
+                if v1 != v2:
+                    logging.warning('Field "%s" has Width=%d, expected %d',
+                        fldName, v1, v2)
+                    ok = False
+
+            if 'Default' in fld:
+                v1 = defn.GetDefault()
+                v2 = fld['Default']
+                if v1 != v2:
+                    logging.warning('Field "%s" has Default="%s", expected "%s"',
+                        fldName, v1, v2)
+                    ok = False
+
+            if 'Nullable' in fld:
+                v1 = bool(defn.IsNullable())
+                v2 = fld['Nullable']
+                if v1 != v2:
+                    logging.warning('Field "%s" has Nullable=%s, expected %s',
+                        fldName, v1, v2)
+                    ok = False
+
+            if 'Unique' in fld:
+                v1 = bool(defn.IsUnique())
+                v2 = fld['Unique']
+                if v1 != v2:
+                    logging.warning('Field "%s" has Unique=%s, expected %s',
+                        fldName, v1, v2)
+                    ok = False
+
+    return ok
+
+# Create output layer
+def createOutputLayer(ds, layername, srs=None, options=None):
+    logging.info('Creating new destination layer "%s"', layername)
+    geom_type = options['geometry-type']
+    lco = options.get('options', None)
+
+    drv = ds.GetDriver()
+    if geom_type != ogr.wkbNone and drv.ShortName == 'PostgreSQL':
+        # “Important to set to 2 for 2D layers as it has constraints on the geometry
+        # dimension during loading.”
+        # — https://gdal.org/drivers/vector/pg.html#layer-creation-options
+        if ogr.GT_HasM(geom_type):
+            if ogr.GT_HasZ(geom_type):
+                dim = 'XYZM'
+            else:
+                dim = 'XYM'
+        elif ogr.GT_HasZ(geom_type):
+            dim = '3'
+        else:
+            dim = '2'
+        if lco is None:
+            lco = []
+        lco = ['dim=' + dim] + lco # prepend dim=
+
+    kwargs = { 'geom_type': geom_type }
+    if srs is not None:
+        kwargs['srs'] = srs
+    if lco is not None:
+        kwargs['options'] = lco
+    logging.debug('CreateLayer(%s, geom_type="%s"%s%s)', layername,
+        ogr.GeometryTypeToName(geom_type),
+        ', srs="' + kwargs['srs'].GetName() + '"' if 'srs' in kwargs else '',
+        ', options=' + str(kwargs['options']) if 'options' in kwargs else '')
+    lyr = dso.CreateLayer(layername, **kwargs)
+    if lyr is None:
+        raise Exception(f'Could not create destination layer "{layername}"')
+
+    fields = options['fields']
+    if len(fields) > 0 and not lyr.TestCapability(ogr.OLCCreateField):
+        raise Exception(f'Destination layer "{layername}" lacks field creation capability')
+
+    # set up output schema
+    for fld in fields:
+        fldName = fld['Name']
+        defn = ogr.FieldDefn()
+        defn.SetName(fldName)
+
+        if 'AlternativeName' in fld:
+            v = fld['AlternativeName']
+            logging.debug('Set AlternativeName=%s on output field "%s"', str(v), fldName)
+            defn.SetAlternativeName(v)
+
+        if 'Comment' in fld:
+            v = fld['Comment']
+            logging.debug('Set Comment=%s on output field "%s"', str(v), fldName)
+            defn.SetComment(v)
+
+        if 'Type' in fld:
+            v = fld['Type']
+            logging.debug('Set Type=%d (%s) on output field "%s"',
+                v, ogr.GetFieldTypeName(v), fldName)
+            defn.SetType(v)
+
+        if 'SubType' in fld:
+            v = fld['SubType']
+            logging.debug('Set SubType=%d (%s) on output field "%s"',
+                v, ogr.GetFieldSubTypeName(v), fldName)
+            defn.SetSubType(v)
+
+        if 'TZFlag' in fld:
+            v, n = fld['TZFlag']
+            logging.debug('Set TZFlag=%d (%s) on output field "%s"', v, n, fldName)
+            defn.SetTZFlag(v)
+
+        if 'Precision' in fld:
+            v = fld['Precision']
+            logging.debug('Set Precision=%d on output field "%s"', v, fldName)
+            defn.SetPrecision(v)
+
+        if 'Width' in fld:
+            v = fld['Width']
+            logging.debug('Set Width=%d on output field "%s"', v, fldName)
+            defn.SetWidth(v)
+
+        if 'Default' in fld:
+            v = fld['Default']
+            logging.debug('Set Default=%s on output field "%s"', v, fldName)
+            defn.SetDefault(v)
+
+        if 'Nullable' in fld:
+            v = fld['Nullable']
+            logging.debug('Set Nullable=%s on output field "%s"', v, fldName)
+            defn.SetNullable(v)
+
+        if 'Unique' in fld:
+            v = fld['Unique']
+            logging.debug('Set Unique=%s on output field "%s"', v, fldName)
+            defn.SetUnique(v)
+
+        if lyr.CreateField(defn, approx_ok=False) != GDAL_CE_None:
+            raise Exception('Could not create field "{fldName}"')
+        logging.debug('Added field "%s" to output layer "%s"', fldName, layername)
+
+    # flush before calling StartTransaction() so we're not tryingn to
+    # rollback changes on a non-existing table
+    lyr.SyncToDisk()
+    return lyr
+
+# Setup output field mapping, modifying the sources dictionary in place.
+def setOutputFieldMap(defn, sources):
+    fieldMap = {}
+    n = defn.GetFieldCount()
+    for i in range(n):
+        fld = defn.GetFieldDefn(i)
+        fldName = fld.GetName()
+        fieldMap[fldName] = i
+
+    for source in sources:
+        src = source['source']['path']
+        fieldMap2 = source['import'].get('fields', None)
+        if fieldMap2 is None:
+            fieldMap2 = {}
+        else:
+            if isinstance(fieldMap2, list):
+                # convert list to identity dictionary
+                fieldMap2 = { fld: fld for fld in fieldMap2 }
+
+            for ifld, ofld in fieldMap2.items():
+                i = fieldMap.get(ofld, None)
+                if i is None:
+                    raise Exception(f'Ouput layer has no field named "{ofld}"')
+                fieldMap2[ifld] = i
+
+        source['import']['fields'] = fieldMap2
+
+    return fieldMap
+
+# 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('"', '""') + '"'
+
+# Clear the given layer (wipe all its features)
+def clearLayer(ds, lyr):
+    n = -1
+    if lyr.TestCapability(ogr.OLCFastFeatureCount):
+        n = lyr.GetFeatureCount(force=0)
+        if n == 0:
+            # nothing to clear, we're good
+            return
+    layername_esc = escapeIdentifier(lyr.GetName())
+
+    # GDAL <3.9 doesn't have lyr.GetDataset() so we pass the DS along with the layer
+    drv = ds.GetDriver()
+    if drv.ShortName == 'PostgreSQL':
+        # https://www.postgresql.org/docs/15/sql-truncate.html
+        query = 'TRUNCATE TABLE {table} CONTINUE IDENTITY RESTRICT'
+        op = 'Truncating'
+    else:
+        query = 'DELETE FROM {table}'
+        op = 'Clearing'
+    logging.info('%s table %s (former feature count: %s)', op,
+        layername_esc, str(n) if n >= 0 else 'unknown')
+    ds.ExecuteSQL(query.format(table=layername_esc))
+
+# Extract an archive file into the given destination directory.
+def extractArchive(path, destdir, fmt=None, patterns=None, exact_matches=None):
+    if fmt is None:
+        suffix = path.suffix
+        if suffix is None or suffix == '' or not suffix.startswith('.'):
+            raise Exception(f'Could not infer archive format from "{path}"')
+        fmt = suffix.removeprefix('.')
+
+    fmt = fmt.lower()
+    logging.debug('Unpacking %s archive %s into %s', fmt, path, destdir)
+
+    if fmt == 'zip':
+        from zipfile import ZipFile
+        logging.debug('Opening %s as ZipFile', path)
+        with ZipFile(path, mode='r') as z:
+            namelist = listArchiveMembers(z.namelist(),
+                patterns=patterns, exact_matches=exact_matches)
+            z.extractall(path=destdir, members=namelist)
+    else:
+        raise Exception(f'Unknown archive format "{fmt}"')
+
+# List archive members matching the given parterns and/or exact matches
+def listArchiveMembers(namelist, patterns, exact_matches=None):
+    if patterns is None and exact_matches is None:
+        # if neither patterns nor exact_matches are given we'll extract the entire archive
+        return namelist
+    if patterns is None:
+        patterns = []
+    if exact_matches is None:
+        exact_matches = []
+
+    members = []
+    for name in namelist:
+        ok = False
+        if name in exact_matches:
+            # try exact matches first
+            logging.debug('Listed archive member %s (exact match)', name)
+            members.append(name)
+            ok = True
+            continue
+        # if there are no exact matches, try patterns one by one in the supplied order
+        for pat in patterns:
+            if fnmatchcase(name, pat):
+                logging.debug('Listed archive member %s (matching pattern "%s")', name, pat)
+                members.append(name)
+                ok = True
+                break
+        if not ok:
+            logging.debug('Ignoring archive member %s', name)
+    return members
+
+# Import a source layer
+def importSource(lyr, path=None, unar=None, args={}, cachedir=None, extent=None):
+    if unar is None:
+        return importSource2(lyr, str(path), args=args,
+                    basedir=cachedir, extent=extent)
+
+    if cachedir is not None:
+        path = cachedir.joinpath(path)
+
+    ds_srcpath = Path(args['path'])
+    if ds_srcpath.is_absolute():
+        # treat absolute paths as relative to the archive root
+        logging.warning('%s is absolute, removing leading anchor', ds_srcpath)
+        ds_srcpath = ds_srcpath.relative_to(ds_srcpath.anchor)
+    ds_srcpath = str(ds_srcpath)
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        logging.debug('Created temporary directory %s', tmpdir)
+        extractArchive(path, tmpdir,
+            fmt=unar.get('format', None),
+            patterns=unar.get('patterns', None),
+            exact_matches=[ds_srcpath])
+        return importSource2(lyr, ds_srcpath, args=args,
+                    basedir=Path(tmpdir), extent=extent)
+
+# Import a source layer (already extracted)
+# This is more or less like ogr2ogr/GDALVectorTranslate() but we roll
+# out our own (slower) version because GDALVectorTranslate() insists in
+# calling StartTransaction() https://github.com/OSGeo/gdal/issues/3403
+# while we want a single transaction for the entire desination layer,
+# including truncation, source imports, and metadata changes.
+def importSource2(lyr_dst, path, args={}, basedir=None, extent=None):
+    kwargs, _ = setOpenExArgs(args, flags=GDAL_OF_READONLY|GDAL_OF_VERBOSE_ERROR)
+    path2 = path if basedir is None else str(basedir.joinpath(path))
+
+    logging.debug('OpenEx(%s, %s)', path2, str(kwargs))
+    ds = gdal.OpenEx(path2, **kwargs)
+    if ds is None:
+        raise Exception(f'Could not open {path2}')
+
+    layername = args.get('layername', None)
+    if layername is None:
+        idx = 0
+        lyr = ds.GetLayerByIndex(idx)
+        msg = '#' + str(idx)
+        if lyr is not None:
+            layername = lyr.GetName()
+            msg += ' ("' + layername + '")'
+    else:
+        lyr = ds.GetLayerByName(layername)
+        msg = '"' + layername + '"'
+    if lyr is None:
+        raise Exception(f'Could not get requested layer {msg} from {path2}')
+
+    logging.info('Importing layer %s from "%s"', msg, path)
+    canIgnoreFields = lyr.TestCapability(ogr.OLCIgnoreFields)
+
+    defn = lyr.GetLayerDefn()
+    n = defn.GetGeomFieldCount()
+    if n != 1: # XXX
+        logging.warning('Source layer "%s" has %d != 1 geometry fields', layername, n)
+
+    n = defn.GetFieldCount()
+    fieldMap = [-1] * n
+    fields = args['fields']
+    fieldSet = set()
+
+    for i in range(n):
+        fld = defn.GetFieldDefn(i)
+        fldName = fld.GetName()
+        fieldMap[i] = v = fields.get(fldName, -1)
+        fieldSet.add(fldName)
+
+        if v < 0 and canIgnoreFields:
+            # call SetIgnored() on unwanted source fields
+            logging.debug('Set Ignored=True on output field "%s"', fldName)
+            fld.SetIgnored(True)
+    defn = None
+
+    logging.debug('Field map: %s', str(fieldMap))
+    for fld in fields:
+        if not fld in fieldSet:
+            logging.warning('Source layer "%s" has no field named "%s", ignoring', layername, fld)
+
+    count0 = -1
+    if lyr.TestCapability(ogr.OLCFastFeatureCount):
+        count0 = lyr.GetFeatureCount(force=0)
+
+    count1 = -1
+    if args.get('spatial-filter', True) and extent is not None:
+        # TODO reprojection
+        lyr.SetSpatialFilterRect(*extent)
+        if lyr.TestCapability(ogr.OLCFastFeatureCount):
+            count1 = lyr.GetFeatureCount(force=0)
+
+    if count0 >= 0:
+        if count1 >= 0:
+            logging.info('Source layer "%s" has %d features (of which %d within extent)',
+                layername, count0, count1)
+        else:
+            logging.info('Source layer "%s" has %d features', layername, count0)
+
+    #print(lyr.GetSupportedSRSList())
+    # TODO reprojection
+
+    defn_dst = lyr_dst.GetLayerDefn()
+    eGType_dst = defn_dst.GetGeomType()
+
+    n = 0
+    #mismatch = 0
+    feature = lyr.GetNextFeature()
+    while feature is not None:
+        #print(lyr.GetFeaturesRead())
+        feature2 = ogr.Feature(defn_dst)
+        # TODO MakeValid, force geometry, dimension
+        feature2.SetFromWithMap(feature, False, fieldMap)
+
+        geom = feature2.GetGeometryRef()
+        eGType = geom.GetGeometryType()
+        if eGType != eGType_dst:
+            if eGType == ogr.wkbPolygon and eGType_dst == ogr.wkbMultiPolygon:
+                geom = ogr.ForceToMultiPolygon(geom)
+            else:
+                raise Exception('Not implemented')
+            feature2.SetGeometryDirectly(geom)
+            #mismatch += 1
+
+        lyr_dst.CreateFeature(feature2)
+        n += 1
+        feature = lyr.GetNextFeature()
+
+    #print(mismatch)
+    lyr = None
+    logging.info('%d features imported from source layer "%s"', n, layername)
+
+
+if __name__ == '__main__':
+    logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO)
+
+    parser = argparse.ArgumentParser(description='Extract and import GIS layers.')
+    parser.add_argument('--cachedir', default=None,
+        help=f'cache directory (default: {os.curdir})')
+    parser.add_argument('--debug', action='count', default=0,
+        help=argparse.SUPPRESS)
+    parser.add_argument('groupname', nargs='*', help='group layer name(s) to process')
+    args = parser.parse_args()
+
+    if args.debug > 0:
+        logging.getLogger().setLevel(logging.DEBUG)
+    if args.debug > 1:
+        gdal.ConfigurePythonLogging(enable_debug=True)
+
+    common.load_config(groupnames=None if args.groupname == [] else args.groupname)
+
+    # validate configuration
+    if 'dataset' not in common.config:
+        raise Exception('Configuration does not specify output dataset')
+
+    layers = common.config.get('layers', {})
+    for layername, layerdefs in layers.items():
+        for idx, layerdef in enumerate(layerdefs['sources']):
+            importdef = layerdef.get('import', None)
+            if importdef is None:
+                raise Exception(f'Output layer "{layername}" source #{idx} has no import definition')
+
+            sourcedef = layerdef.get('source', None)
+            unar = None if sourcedef is None else sourcedef.get('unar', None)
+            src = None if sourcedef is None else sourcedef['cache'].get('path', None)
+
+            ds_srcpath = importdef.get('path', None)
+            if src is None and unar is None and ds_srcpath is not None:
+                # fallback to importe:path if there is no unarchiving receipe
+                src = Path(ds_srcpath)
+            if unar is not None and ds_srcpath is None:
+                raise Exception(f'Output layer "{layername}" source #{idx} has no import source path')
+            if src is None:
+                raise Exception(f'Output layer "{layername}" source #{idx} has no source path')
+            layerdef['source'] = { 'path': src, 'unar': unar }
+
+    # set global GDAL/OGR configuration options
+    for pszKey, pszValue in common.config.get('GDALconfig', {}).items():
+        logging.debug('gdal.SetConfigOption(%s, %s)', pszKey, pszValue)
+        gdal.SetConfigOption(pszKey, pszValue)
+
+    # open output dataset (possibly create it first)
+    dso = openOutputDS(common.config['dataset'])
+
+    validateSchema(layers,
+        drvo=dso.GetDriver(),
+        lco_defaults=common.config['dataset'].get('create-layer-options', None))
+
+    # get configured Spatial Reference System and extent
+    srs = getSRS(common.config.get('SRS', None))
+    extent = common.config.get('extent', None)
+    if extent is not None:
+        logging.debug('Configured extent in %s: %s',
+            srs.GetName(), ', '.join(map(str, extent)))
+
+    cachedir = Path(args.cachedir) if args.cachedir is not None else None
+    rv = 0
+    for layername, layerdef in layers.items():
+        sources = layerdef['sources']
+        if sources is None or len(sources) < 1:
+            logging.warning('Output layer "%s" has no definition, skipping', layername)
+            continue
+
+        logging.info('Processing output layer "%s"', layername)
+        transaction = False
+        try:
+            # get output layer
+            outLayerIsNotEmpty = True
+            lco = layerdef.get('create', None)
+            lyr = dso.GetLayerByName(layername)
+            if lyr is not None:
+                # TODO dso.DeleteLayer(layername) if --overwrite and dso.TestCapability(ogr.ODsCDeleteLayer)
+                # (Sets OVERWRITE=YES for PostgreSQL and GPKG.)
+                validateOutputLayer(lyr, srs=srs, options=lco)
+                # TODO bail out if all source files are older than lyr's last_change
+            elif not dso.TestCapability(ogr.ODsCCreateLayer):
+                raise Exception(f'Output driver {dso.GetDriver().ShortName} does not support layer creation')
+            elif lco is None or len(lco) < 1:
+                raise Exception(f'Missing schema for new output layer "{layername}"')
+            else:
+                lyr = createOutputLayer(dso, layername, srs=srs, options=lco)
+                outLayerIsNotEmpty = False
+
+            if not lyr.TestCapability(ogr.OLCSequentialWrite):
+                raise Exception(f'Output layer "{layername}" has no working CreateFeature() method')
+
+            # setup output field mapping in the sources dictionary
+            setOutputFieldMap(lyr.GetLayerDefn(), sources)
+
+            # start transaction if possible
+            if lyr.TestCapability(ogr.OLCTransactions):
+                logging.debug('Starting transaction')
+                transaction = lyr.StartTransaction() == ogr.OGRERR_NONE
+            else:
+                logging.warning('Unsafe update, output layer "%s" does not support transactions', layername)
+
+            if outLayerIsNotEmpty:
+                # clear non-empty output layer
+                clearLayer(dso, lyr)
+
+            if lyr.SetMetadataItem('DESCRIPTION', layerdef.get('description', '')) != GDAL_CE_None:
+                logging.warning('Could not set description metadata')
+
+            for source in sources:
+                # import source layers
+                importSource(lyr, **source['source'], args=source['import'],
+                    cachedir=cachedir, extent=extent)
+
+            if transaction:
+                # commit transaction
+                logging.debug('Committing transaction')
+                transaction = False
+                if lyr.CommitTransaction() != ogr.OGRERR_NONE:
+                    logging.error('Could not commit transaction')
+                    rv = 1
+
+        except Exception:
+            if transaction:
+                logging.error('Exception occured in transaction, rolling back')
+                try:
+                    if lyr.RollbackTransaction() != ogr.OGRERR_NONE:
+                        logging.error('Could not rollback transaction')
+                except RuntimeError:
+                    logging.exception('Could not rollback transaction')
+            logging.exception('Could not import layer "%s"', layername)
+            rv = 1
+
+        finally:
+            # close output layer
+            lyr = None
+
+    dso = None
+    srs = None
+    exit(rv)