#!/usr/bin/python3 #---------------------------------------------------------------------- # Backend utilities for the Klimatanalys Norr project (extract/import layers) # Copyright © 2024 Guilhem Moulin # # 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 . #---------------------------------------------------------------------- import os import logging import argparse import tempfile import re import math 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 parseFieldType(name): if name is None: raise Exception('parseFieldType(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 parseSubFieldType(name): if name is None: raise Exception('parseSubFieldType(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 "{name}"') # Parse timezone TZ_RE = re.compile(r'(?:UTC\b)?([\+\-]?)([0-9][0-9]):?([0-9][0-9])', flags=re.IGNORECASE) def parseTimeZone(tz): if tz is None: raise Exception('parseTimeZone(None)') tz2 = tz.lower() if tz2 == 'none': return ogr.TZFLAG_UNKNOWN elif tz2 == 'local': return ogr.TZFLAG_LOCALTIME elif tz2 == 'utc' or tz2 == 'gmt': return ogr.TZFLAG_UTC m = TZ_RE.fullmatch(tz) if m is None: raise Exception(f'Invalid timezone "{tz}"') tzSign = m.group(1) tzHour = int(m.group(2)) tzMinute = int(m.group(3)) if tzHour > 14 or tzMinute >= 60 or tzMinute % 15 != 0: raise Exception(f'Invalid timezone "{tz}"') tzFlag = tzHour*4 + int(tzMinute/15) if tzSign == '-': tzFlag = 100 - tzFlag else: tzFlag += 100 return tzFlag # Pretty-print timezone flag, cf. # ogr/ogrutils.cpp:OGRGetISO8601DateTime() def formatTZFlag(tzFlag): if tzFlag is None: raise Exception('printTimeZone(None)') if tzFlag == ogr.TZFLAG_UNKNOWN: return 'none' elif tzFlag == ogr.TZFLAG_LOCALTIME: return 'local' elif tzFlag == ogr.TZFLAG_UTC: return 'UTC' tzOffset = abs(tzFlag - 100) * 15; tzHour = int(tzOffset / 60); tzMinute = int(tzOffset % 60); tzSign = '+' if tzFlag > 100 else '-' return f'{tzSign}{tzHour:02}{tzMinute:02}' # 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': # (WARN support added in GDAL 3.7) fld_def2['Comment'] = v elif k2 == 'type': fld_def2['Type'] = parseFieldType(v) elif k2 == 'subtype': fld_def2['SubType'] = parseSubFieldType(v) elif k2 == 'tz': # (WARN support added in GDAL 3.8) fld_def2['TZFlag'] = parseTimeZone(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 # Convert extent [minX, minY, maxX, maxY] into a polygon and assign the # given SRS. Like apps/ogr2ogr_lib.cpp, we segmentize 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 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))) 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) polygon.TransformTo(srs) segment_distance_metre = 10 * 1000 if srs.IsGeographic(): dfMaxLength = segment_distance_metre / math.radians(srs.GetSemiMajor()) polygon.Segmentize(dfMaxLength) elif srs.IsProjected(): dfMaxLength = segment_distance_metre / srs.GetLinearUnits() polygon.Segmentize(dfMaxLength) return polygon # 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 = fld['TZFlag'] if v1 != v2: logging.warning('Field "%s" has TZFlag=%d (%s), expected %d (%s)', fldName, v1, formatTZFlag(v1), v2, formatTZFlag(v2)) 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 = fld['TZFlag'] logging.debug('Set TZFlag=%d (%s) on output field "%s"', v, formatTZFlag(v), 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()) # XXX 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) srs = lyr.GetSpatialRef() if srs is None: raise Exception('Source layer has no SRS') srs_dst = lyr_dst.GetSpatialRef() if srs_dst.IsSame(srs): logging.debug('Both source and destination have the same SRS (%s), skipping coordinate transformation', srs_dst.GetName()) ct = None else: # TODO Use GetSupportedSRSList() and SetActiveSRS() with GDAL ≥3.7.0 # when possible, see apps/ogr2ogr_lib.cpp logging.debug('Creating transforming from source SRS (%s) to destination SRS (%s)', srs.GetName(), srs_dst.GetName()) ct = osr.CoordinateTransformation(srs, srs_dst) if ct is None: raise Exception(f'Could not create transformation from source SRS ({srs.GetName()}) ' + f'to destination SRS ({srs_dst.GetName()})') defn = lyr.GetLayerDefn() n = defn.GetGeomFieldCount() if n != 1: # TODO Add support for multiple geometry fields 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: if extent.GetSpatialReference().IsSame(srs): extent2 = extent else: extent2 = extent.Clone() if extent2.TransformTo(srs) != ogr.OGRERR_NONE: raise Exception(f'Could not transform extent {extent.ExportToWkt()} to {srs.GetName()}') #logging.debug('Applying extent: %s', extent2.ExportToWkt()) lyr.SetSpatialFilter(extent2) 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) defn_dst = lyr_dst.GetLayerDefn() eGType_dst = defn_dst.GetGeomType() eGType_dst_HasZ = ogr.GT_HasZ(eGType_dst) eGType_dst_HasM = ogr.GT_HasM(eGType_dst) dGeomIsUnknown = ogr.GT_Flatten(eGType_dst) == ogr.wkbUnknown n = 0 mismatch = {} feature = lyr.GetNextFeature() while feature is not None: feature2 = ogr.Feature(defn_dst) feature2.SetFromWithMap(feature, False, fieldMap) geom = feature2.GetGeometryRef() if ct is not None and geom.Transform(ct) != ogr.OGRERR_NONE: raise Exception('Could not apply coordinate transformation') eGType = geom.GetGeometryType() if eGType != eGType_dst and not dGeomIsUnknown: # Promote to multi, cf. apps/ogr2ogr_lib.cpp:ConvertType() eGType2 = eGType if eGType == ogr.wkbTriangle or eGType == ogr.wkbTIN or eGType == ogr.wkbPolyhedralSurface: eGType2 = ogr.wkbMultiPolygon elif not ogr.GT_IsSubClassOf(eGType, ogr.wkbGeometryCollection): eGType2 = ogr.GT_GetCollection(eGType) eGType2 = ogr.GT_SetModifier(eGType2, eGType_dst_HasZ, eGType_dst_HasM) if eGType2 == eGType_dst: mismatch[eGType] = mismatch.get(eGType, 0) + 1 geom = ogr.ForceTo(geom, eGType_dst) # TODO call MakeValid()? else: raise Exception(f'Conversion from {ogr.GeometryTypeToName(eGType)} ' f'to {ogr.GeometryTypeToName(eGType_dst)} not implemented') feature2.SetGeometryDirectly(geom) if lyr_dst.CreateFeature(feature2) != ogr.OGRERR_NONE: raise Exception(f'Could not transfer source feature #{feature.GetFID()}') n += 1 feature = lyr.GetNextFeature() lyr = None logging.info('Imported %d features from source layer "%s"', n, layername) if len(mismatch) > 0: mismatches = [ str(n) + '× ' + ogr.GeometryTypeToName(t) for t,n in sorted(mismatch.items(), key=lambda x: x[1]) ] logging.info('Forced conversion to %s: %s', ogr.GeometryTypeToName(eGType_dst), ', '.join(mismatches)) if __name__ == '__main__': common.init_logger(app=os.path.basename(__file__), 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 = getExtent(common.config.get('extent', None), srs=srs) 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) description = layerdef.get('description', None) if description is not None and lyr.SetMetadataItem('DESCRIPTION', 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 extent = None exit(rv)