1#  coding=utf-8
   2#
   3#  Author: Ernesto Arredondo Martinez (ernestone@gmail.com)
   4#  Created: 7/6/19 18:23
   5#  Last modified: 7/6/19 14:27
   6#  Copyright (c) 2019
   7"""
   8.. include:: ../README.md
   9"""
  10
  11import logging
  12import os
  13import re
  14from collections import OrderedDict, namedtuple
  15from logging import Logger
  16
  17import math
  18from osgeo import gdal, ogr, osr
  19from osgeo.ogr import ODsCCreateLayer, OLCAlterFieldDefn, OLCCreateField, ODsCTransactions, \
  20    ODsCDeleteLayer, OLCTransactions, Geometry, ODrCCreateDataSource, GeomFieldDefn
  21
  22from apb_extra_utils import misc as utils, utils_logging
  23
  24DRVR_PG = 'PostgreSQL'
  25DRVR_CSV = "CSV"
  26DRVR_GPKG = "GPKG"
  27DRVR_GEOJSON = "GeoJSON"
  28DRVR_KML = "KML"
  29
  30SIMPLIFY_TOLERANCE = 1e-6
  31DRIVERS_GDAL_NOT_DELETE_LAYERS = (DRVR_PG,)
  32DRIVERS_OLD_SRS_AXIS_MAPPING_4326 = (DRVR_GEOJSON, DRVR_GPKG, DRVR_PG)
  33DRIVERS_NO_SPATIAL_IDX = (DRVR_GEOJSON, DRVR_CSV, DRVR_PG)
  34PREFFIX_GEOMS_LAYERS_GDAL = 'geom_'
  35PREFFIX_GEOMS_LAYERS_GDAL_CSV = '_WKT'
  36
  37# Cache global: {ds_name: {'count': int, 'layers': set[str]}}
  38_CACHE_DS_POSTGIS: dict[str, ogr.DataSource] = {}
  39_CACHE_NOM_LAYERS_DS_GDAL = {}
  40
  41print_debug = logging.debug
  42print_warning = logging.warning
  43
  44
  45def get_geom_field_defn(layer_feature_defn, id_geom: int = 0):
  46    """
  47    Retorna el GeomFieldDefn indicado solo cuando existe al menos un campo geométrico.
  48
  49    Args:
  50        layer_feature_defn (osgeo.ogr.FeatureDefn): una instancia de FeatureDefn (vease https://gdal.org/en/stable/api/python/vector_api.html#featuredefn)
  51        id_geom (int=0): índice del campo geométrico a recuperar.
  52
  53    Returns:
  54        osgeo.ogr.GeomFieldDefn: vease https://gdal.org/en/stable/api/python/vector_api.html#geomfielddefn
  55    """
  56    if not layer_feature_defn:
  57        return None
  58
  59    if id_geom < 0:
  60        raise ValueError(f"id_geom debe ser >= 0. Valor recibido: {id_geom}")
  61
  62    count_geoms = layer_feature_defn.GetGeomFieldCount()
  63    if count_geoms == 0:
  64        return None
  65
  66    if id_geom >= count_geoms:
  67        raise ValueError(
  68            f"id_geom fuera de rango ({id_geom}). Geometrías disponibles: {count_geoms} (índices 0..{count_geoms - 1})"
  69        )
  70
  71    return layer_feature_defn.GetGeomFieldDefn(id_geom)
  72
  73
  74def srs_ref_from_epsg_code(code_epsg: int, old_axis_mapping=False) -> osr.SpatialReference:
  75    """
  76
  77    Args:
  78        code_epsg (int):
  79        old_axis_mapping (bool=False): Por cambio entre GDAL2.0 y GDAL3.0 se (vease
  80    https://github.com/OSGeo/gdal/blob/master/gdal/MIGRATION_GUIDE.TXT) por defecto se utiliza el mapeo de ejes
  81    clasico LONG,LAT
  82
  83    Returns:
  84        srs (osr.SpatialReference)
  85    """
  86    srs = osr.SpatialReference()
  87    if old_axis_mapping:
  88        try:
  89            from osgeo.osr import OAMS_TRADITIONAL_GIS_ORDER
  90            print_debug("OAMS_TRADITIONAL_GIS_ORDER")
  91            srs.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER)
  92        except ImportError:
  93            print_warning("OAMS_TRADITIONAL_GIS_ORDER not available")
  94
  95    ret = srs.ImportFromEPSG(code_epsg)
  96    if ret != 0:
  97        raise Warning("No se puede retornar un osgeo.osr.SpatialReference EPSG para el codigo '{}'!!".format(code_epsg))
  98
  99    return srs
 100
 101
 102def layer_gdal_from_file(path_file: str,
 103                         nom_driver: str = DRVR_GEOJSON,
 104                         nom_geom: str = None,
 105                         default_srs_epsg_code: int = 4326,
 106                         default_order_long_lat: bool = True):
 107    """
 108
 109    Args:
 110        path_file (str):
 111        nom_driver (str=DRVR_GEOJSON):
 112        nom_geom (str=None): si se informa devolverá la layer solo con la geometria especificada
 113        default_srs_epsg_code (int=4326): codigo del sistema de coordenadas que asignará por defecto si la layer
 114                                          NO tiene sistema definido
 115        default_order_long_lat (bool=True): si no viene informado el SRS y se usa el default EPSG4326 entonces se supone
 116                                          por defecto que las geometria vienen en orden antiguo LONGITUD, LATITUD
 117
 118    Returns:
 119        layer_gdal (osgeo.ogr.Layer), nom_layer (str), datasource_gdal (osgeo.ogr.DataSource)
 120    """
 121    drvr = ogr.GetDriverByName(nom_driver)
 122
 123    nom_layer, ext = utils.split_ext_file(os.path.basename(path_file))
 124    if ext.lower().endswith("zip"):
 125        path_file = "/vsizip/{}".format(path_file)
 126
 127    ds_vector_file = drvr.Open(path_file, 0)
 128
 129    a_layer = None
 130    if ds_vector_file:
 131        a_layer = ds_vector_file.GetLayer(0)
 132
 133        if nom_geom:
 134            lay_def = a_layer.GetLayerDefn()
 135            nom_geom = nom_geom.strip().upper()
 136
 137            if lay_def.GetGeomFieldCount() > 1:
 138                idx_geom = lay_def.GetGeomFieldIndex(nom_geom)
 139                if idx_geom < 0:
 140                    raise Exception("El fichero vectorial '{}' no contiene una geometria con el nombre '{}'".format(
 141                        ds_vector_file.GetName(), nom_geom))
 142                elif idx_geom > 0:
 143                    # Convertimos a layer en MEMORY para poder cambiar estructura
 144                    ds_mem = ds_gdal_memory()
 145                    a_layer = create_layer_from_layer_gdal_on_ds_gdal(ds_mem, a_layer, nom_layer, nom_geom,
 146                                                                      exclude_cols_geoms=False, null_geoms=True)
 147                    if a_layer:
 148                        ds_vector_file = ds_mem
 149
 150            elif not a_layer.GetGeometryColumn() and (geom_def := get_geom_field_defn(lay_def)):
 151                geom_def.SetName(nom_geom)
 152
 153        if a_layer:
 154            # Por si no carga el SRS se asigna el :default_srs_epsg_code (defecto epsg:4326)
 155            for gf in geom_fields_layer_gdal(a_layer):
 156                if not gf.GetSpatialRef():
 157                    gf.SetSpatialRef(srs_ref_from_epsg_code(default_srs_epsg_code, default_order_long_lat))
 158
 159    return a_layer, nom_layer, ds_vector_file
 160
 161
 162def datasource_gdal_vector_file(nom_driver_gdal, nom_ds, a_dir, create=None, from_zip=False, **create_options):
 163    """
 164    Crea datasource gdal para driver de tipo fichero
 165    Args:
 166        nom_driver_gdal (str):
 167        nom_ds (str):
 168        a_dir (str):
 169        create (bool=None): Por defecto crea el datasource si este no existe y se abre sin problemas previamente.
 170                        Si False entonces NO lo crea en ningún caso, y si True lo creará sin intentar abrirlo
 171        from_zip (bool=False):
 172        **create_options: lista claves valores con opciones de creacion para el datasource de GDAL
 173                        p.e.: NameField="SEQENTITAT", DescriptionField="SW_URN"
 174
 175    Returns:
 176        datasource_gdal (osgeo.ogr.DataSource), overwrited (bool)
 177    """
 178    driver, exts_driver = driver_gdal(nom_driver_gdal)
 179    if not exts_driver:
 180        raise Exception("ERROR! - El driver GDAL {} no es de tipo fichero vectorial".format(driver))
 181
 182    base_path_file = os.path.normpath(os.path.join(a_dir, nom_ds.strip().lower()))
 183    _, ext_base = os.path.splitext(base_path_file)
 184
 185    vsi_prefix = ""
 186    if from_zip:
 187        ext_file = "zip"
 188        vsi_prefix = "/vsizip/"
 189    else:
 190        if "topojson" in (ext.lower() for ext in exts_driver):
 191            ext_file = "topojson"
 192        else:
 193            ext_file = exts_driver[0]
 194
 195    if (os.path.exists(base_path_file) and os.path.isfile(
 196            base_path_file)) or ext_base.lower() == f'.{ext_file.lower()}':
 197        path_file = base_path_file
 198    else:
 199        path_file = "{}.{}".format(base_path_file, ext_file)
 200
 201    overwrited = False
 202    datasource_gdal = None
 203    if not create and os.path.exists(path_file):
 204        open_path_file = "{}{}".format(vsi_prefix, path_file)
 205        datasource_gdal = driver.Open(open_path_file, 1)
 206        if datasource_gdal is None:
 207            datasource_gdal = driver.Open(open_path_file)
 208
 209    if datasource_gdal:
 210        overwrited = True
 211
 212    if create or (create is not False and datasource_gdal is None and driver.TestCapability(ODrCCreateDataSource)):
 213        datasource_gdal = driver.CreateDataSource(path_file,
 214                                                  list("{}={}".format(opt, val) for opt, val in create_options.items()))
 215
 216    return datasource_gdal, overwrited
 217
 218
 219def driver_gdal(nom_driver):
 220    """
 221    Devuelve el driver GDAL solicitado y si es de tipo fichero vectorial tambien devuelve la extension
 222    Si no coincide exactamente devuelve el que tiene nombre más parecido.
 223    Verificar con drivers_ogr_gdal_disponibles() los disponibles
 224
 225    Args:
 226        nom_driver:
 227    Returns:
 228        driver_gdal (osgeo.ogr.Driver), exts_driver (list)
 229    """
 230    driver_gdal = ogr.GetDriverByName(nom_driver)
 231    if driver_gdal is None:
 232        raise ValueError(f"Driver '{nom_driver}' is not available")
 233
 234    metadata = driver_gdal.GetMetadata()
 235    exts_drvr = metadata.get('DMD_EXTENSIONS', "").split(" ") if metadata else []
 236
 237    return driver_gdal, exts_drvr
 238
 239
 240def ds_gdal_memory():
 241    """
 242    Devuelve un osgeo.ogr.DataSource de tipo "memory"
 243
 244    Returns:
 245        datasource_gdal (osgeo.ogr.DataSource)
 246    """
 247    return ogr.GetDriverByName("memory").CreateDataSource("")
 248
 249
 250def set_create_option_list_for_layer_gdal(layer_gdal, drvr_name=DRVR_GPKG, has_geom=True, **extra_opt_list):
 251    """
 252    Devuelve la lista de create options GDAL a partir de una layer, el driver para el que se quiere crear y options
 253    pasadas por defecto
 254
 255    Args:
 256        layer_gdal (ogr.Layer):
 257        drvr_name (str=DRVR_GPKG): nombre de driver GDAL
 258        has_geom (bool=True): Indica si la layer tiene geometria. Si es False no añadirá options de geom como GEOMETRY_NAME, SPATIAL_INDEX, ...
 259        extra_opt_list: lista pares claves-valores que se corresponden con option create list del driver GDAL
 260
 261    Returns:
 262        opt_list (dict)
 263    """
 264    opt_list = set_create_option_list_for_driver_gdal(drvr_name, has_geom=has_geom, **extra_opt_list)
 265
 266    opt_geom_name = "GEOMETRY_NAME"
 267    if has_geom and opt_geom_name not in opt_list:
 268        nom_geom_src = layer_gdal.GetGeometryColumn()
 269        if nom_geom_src:
 270            opt_list[opt_geom_name] = "{}={}".format(opt_geom_name, nom_geom_src)
 271
 272    return opt_list
 273
 274
 275def set_create_option_list_for_driver_gdal(drvr_name=DRVR_GPKG, has_geom=True, **extra_opt_list):
 276    """
 277    Devuelve la lista de create options GDAL a partir de una layer, el driver para el que se quiere crear y options
 278    pasadas por defecto
 279
 280    Args:
 281        drvr_name (str=DRVR_GPKG): nombre de driver GDAL
 282        has_geom (bool=True): Indica si la layer tiene geometria. Si es False no añadirá options de geom como SPATIAL_INDEX, ...
 283        extra_opt_list: lista pares claves-valores que se corresponden con option create list del driver GDAL
 284
 285    Returns:
 286        opt_list (dict)
 287    """
 288    # Se quitan las options que no son de creacion para el driver especificado
 289    drvr, exts_drvr = driver_gdal(drvr_name)
 290    if not drvr:
 291        Exception("!ERROR! - El nombre de driver '{}' no es un driver GDAL disponible".format(drvr_name))
 292
 293    opt_list = {k.upper(): v.upper() for k, v in extra_opt_list.items()}
 294
 295    if not "FID" in opt_list and drvr.GetName() == DRVR_GPKG:
 296        opt_list["FID"] = 'FID=FID_GPKG'
 297
 298    if drvr.GetName() == DRVR_GEOJSON:
 299        if "WRITE_BBOX" not in opt_list:
 300            opt_list["WRITE_BBOX"] = 'WRITE_BBOX=YES'
 301
 302    if drvr.GetName() == DRVR_CSV:
 303        if "CREATE_CSVT" not in opt_list:
 304            opt_list["CREATE_CSVT"] = 'CREATE_CSVT=YES'
 305        if "GEOMETRY" not in opt_list and has_geom:
 306            opt_list["GEOMETRY"] = 'GEOMETRY=AS_WKT'
 307
 308    if drvr.GetMetadataItem('DS_LAYER_CREATIONOPTIONLIST'):
 309        list_opts_drvr = drvr.GetMetadataItem('DS_LAYER_CREATIONOPTIONLIST')
 310        keys_opt_list = list(opt_list.keys())
 311        for n_opt in keys_opt_list:
 312            if list_opts_drvr.find(n_opt) < 0:
 313                opt_list.pop(n_opt)
 314
 315        if has_geom and "SPATIAL_INDEX" not in opt_list and \
 316                list_opts_drvr.find("SPATIAL_INDEX") >= 0 and drvr.GetName() != DRVR_PG:
 317            opt_list["SPATIAL_INDEX"] = 'SPATIAL_INDEX=YES'
 318
 319    return opt_list
 320
 321
 322def copy_layer_gdal_to_ds_gdal(layer_src, ds_gdal_dest, nom_layer=None, nom_geom=None,
 323                               overwrite=True, **extra_opt_list):
 324    """
 325    Copia una layer_gdal a otro datasource gdal
 326
 327    Args:
 328        layer_src:
 329        ds_gdal_dest:
 330        nom_layer (str=None): OPC - Si no viene informado cogerá el nombre de la layer
 331        nom_geom (str=None): OPC - Si se informa se cogerá como el nombre de la geometria de la nueva layer
 332        overwrite (bool=True):
 333        **extra_opt_list (str): Lista claves-valores de opciones para copylayer o createLayer
 334                                del driver de ds_gdal indicado
 335
 336    Returns:
 337        layer_dest (ogr.layer)
 338    """
 339    if not nom_layer:
 340        nom_layer = layer_src.GetName()
 341    nom_layer = nom_layer.strip().lower()
 342
 343    layer_dest = ds_gdal_dest.GetLayerByName(nom_layer)
 344    if layer_dest:
 345        if not overwrite:
 346            return layer_dest
 347        else:
 348            ds_gdal_dest.DeleteLayer(nom_layer)
 349
 350    drvr_name = ds_gdal_dest.GetDriver().GetName()
 351
 352    extra_opt_list["IDENTIFIER"] = "IDENTIFIER={}".format(nom_layer)
 353
 354    lyr_src_has_geom = bool(layer_src.GetGeometryColumn() or nom_geom)
 355
 356    opt_list = set_create_option_list_for_layer_gdal(
 357        layer_src, drvr_name=drvr_name, has_geom=lyr_src_has_geom,
 358        **{str(k).upper(): str(v).upper() for k, v in extra_opt_list.items()})
 359
 360    layer_dest = ds_gdal_dest.CopyLayer(layer_src, nom_layer, list(opt_list.values()))
 361
 362    if layer_dest:
 363        layer_dest_def = layer_dest.GetLayerDefn()
 364        if lyr_src_has_geom and not layer_dest.GetGeometryColumn() and get_geom_field_defn(layer_dest_def):
 365            if not nom_geom:
 366                nom_geom = layer_src.GetGeometryColumn()
 367            else:
 368                nom_geom = nom_geom.upper()
 369
 370            get_geom_field_defn(layer_dest_def).SetName(nom_geom)
 371
 372        if drvr_name == DRVR_GPKG and lyr_src_has_geom:
 373            create_spatial_index_layer_gpkg(ds_gdal_dest, nom_layer)
 374
 375    return layer_dest
 376
 377
 378def layer_gtype_from_geoms(layer_gdal, nom_geom=None):
 379    """
 380    A partir de la 1º geometria informada de una layer_gdal devuelve el tipo de geometria que es. Si no encuentra
 381    devuelve el geom_type de la layer (layer_gdal.GetGeomType())
 382
 383    Args:
 384        layer_gdal (osgeo.ogr.Layer):
 385        nom_geom (str=None): Nombre geometria
 386
 387    Returns:
 388        geom_type (int=layer_gdal.GetGeomType()): Si no encuentra devuelve 0 por defecto (GEOMETRY)
 389    """
 390    idx_geom = 0
 391    if nom_geom:
 392        idx_geom = layer_gdal.GetLayerDefn().GetGeomFieldIndex(nom_geom)
 393
 394    if idx_geom >= 0:
 395        layer_gdal.ResetReading()
 396        return next((f.GetGeomFieldRef(idx_geom).GetGeometryType()
 397                     for f in layer_gdal if f.GetGeomFieldRef(idx_geom)),
 398                    layer_gdal.GetGeomType())
 399
 400
 401def srs_for_layer(layer_src, nom_geom_src=None):
 402    """
 403    Devuelve el SpatialReference de la layer_gdal
 404    Args:
 405        layer_src (ogr.Layer)
 406        nom_geom_src (str=None): Nombre geometria
 407
 408    Returns:
 409        srs_lyr_src (osgeo.osr.SpatialReference)
 410    """
 411    srs_lyr_src = layer_src.GetSpatialRef()
 412    layer_src_def = layer_src.GetLayerDefn()
 413
 414    if act_geom_field := get_geom_field_defn(layer_src_def):
 415        if nom_geom_src:
 416            idx_act_geom_field = layer_src_def.GetGeomFieldIndex(nom_geom_src)
 417            if idx_act_geom_field >= 0:
 418                act_geom_field = layer_src_def.GetGeomFieldDefn(idx_act_geom_field)
 419
 420        if act_geom_field.GetSpatialRef():
 421            srs_lyr_src = act_geom_field.GetSpatialRef()
 422
 423    return srs_lyr_src
 424
 425
 426def create_layer_from_layer_gdal_on_ds_gdal(ds_gdal_dest, layer_src, nom_layer=None, nom_geom_src=None, sel_camps=None,
 427                                            exclude_cols_geoms=True, tolerance_simplify=None, null_geoms=False,
 428                                            gtype_layer_from_geoms=True, epsg_code_dest=None,
 429                                            epsg_code_src_default=4326, old_axis_mapping_srs=None, **extra_opt_list):
 430    """
 431    Crea nuevo layer a partir de layer_gdal.
 432
 433    Args:
 434        ds_gdal_dest (ogr.DataSource):
 435        layer_src (ogr.Layer):
 436        nom_layer (str=None):
 437        nom_geom_src (str=None): Si el nombre no se corresponde con ninguna de las geometrias de la layer_src se utilizará
 438                            como ALIAS de la geometria 0 (la defecto) de la nueva layer resultante
 439        sel_camps (list=None): OPC - Lista de campos a escoger de la layer original
 440        exclude_cols_geoms (bool=True): Por defecto de la lista de columnas alfanuméricas (no geometrias) excluirá las
 441                                       columnas que hagan referencia a alguna de las geometrias
 442        tolerance_simplify (float=None): Tolerancia (distancia minima) en unidades del srs de la layer
 443             Mirar método Simplify() sobre osgeo.ogr.Geometry
 444        null_geoms (bool=False): Por defecto no grabará las filas que la geometria principal (nom_geom) es nula
 445        gtype_layer_from_geoms (bool=True): Por defecto, si gtype de layer origen == 0 (GEOMETRY) deducirá el tipo de
 446                                geometria (POINT, LINE o POLYGON) a partir de la primera informada
 447                                encontrada en la layer_origen
 448        epsg_code_dest (int=None): Codigo EPSG para le que se transformarán las geometrias desde el SRS original
 449        epsg_code_src_default (int=4326): Codigo EPSG que se usará para las layer_src que NO tengan SRS asignado
 450        old_axis_mapping_srs (bool=None): setea si quieres usar old mapping axes en (LONG, LAT) (GDAL >3.2)  
 451        **extra_opt_list (str): Lista claves-valores de opciones para createLayer del driver de ds_gdal indicado
 452
 453    Returns:
 454        layer_out (ogr.Layer)
 455    """
 456    desc_ds_gdal = ds_gdal_dest.GetDescription()
 457    drvr_name = ds_gdal_dest.GetDriver().GetName()
 458    print_debug("Inici crear layer '{}' en ds_gdal '{}'".format(
 459        (nom_layer if nom_layer else layer_src.GetName()).upper(),
 460        desc_ds_gdal if desc_ds_gdal else drvr_name))
 461
 462    geoms_src = geoms_layer_gdal(layer_src)
 463    camps_src = cols_layer_gdal(layer_src)
 464
 465    if nom_geom_src:
 466        nom_geom_src = nom_geom_src.upper()
 467        if len(geoms_src) > 1:
 468            if nom_geom_src not in map(lambda ng: ng.upper(), geoms_src):
 469                raise Exception("Argumento :NOM_GEOM = '{}' erróneo ya que "
 470                                "LAYER_GDAL original no contiene dicha geometria".format(nom_geom_src))
 471        elif len(geoms_src) == 0:
 472            raise Exception("Argumento :NOM_GEOM = '{}' erróneo ya que "
 473                            "LAYER_GDAL original no contiene geometrias".format(nom_geom_src))
 474
 475    if sel_camps:
 476        sel_camps = {ng.upper() for ng in sel_camps}
 477        if not sel_camps.issubset(camps_src):
 478            raise Exception("Argumento :SEL_CAMPS = '[{}]' erróneo ya que "
 479                            "LAYER_GDAL no contiene alguno de los campos indicados".format(",".join(sel_camps)))
 480    else:
 481        sel_camps = set()
 482
 483    if not nom_layer:
 484        nom_layer = layer_src.GetName()
 485    else:
 486        nom_layer = nom_layer.strip()
 487
 488    nom_layer = nom_layer.lower()
 489
 490    gtype = layer_src.GetGeomType()
 491    srs_lyr_src = layer_src.GetSpatialRef()
 492
 493    layer_src_def = layer_src.GetLayerDefn()
 494    act_geom_field = get_geom_field_defn(layer_src_def)
 495    nom_geom_out = None
 496
 497    if act_geom_field:
 498        if nom_geom_src:
 499            idx_act_geom_field = layer_src_def.GetGeomFieldIndex(nom_geom_src)
 500            if idx_act_geom_field >= 0:
 501                act_geom_field = layer_src_def.GetGeomFieldDefn(idx_act_geom_field)
 502                nom_geom_out = fix_affix_geom_name_layer_gdal(nom_geom_src, layer_src)
 503
 504        gtype = act_geom_field.GetType()
 505        if gtype_layer_from_geoms and not gtype:
 506            gtype = layer_gtype_from_geoms(layer_src, nom_geom_src)
 507
 508        if act_geom_field.GetSpatialRef():
 509            srs_lyr_src = act_geom_field.GetSpatialRef()
 510
 511    if sel_camps:
 512        sel_camps = {c.strip().upper() for c in sel_camps}
 513
 514    geoms_src_minus_affix = geoms_layer_gdal(layer_src, PREFFIX_GEOMS_LAYERS_GDAL)
 515    if exclude_cols_geoms:
 516        if sel_camps:
 517            sel_camps.difference_update(geoms_src_minus_affix)
 518        else:
 519            sel_camps = cols_layer_gdal(layer_src).difference(geoms_src_minus_affix)
 520
 521    srs_lyr_dest = None
 522    if not srs_lyr_src and epsg_code_src_default:
 523        srs_lyr_src = srs_ref_from_epsg_code(epsg_code_src_default, old_axis_mapping=old_axis_mapping_srs)
 524    if srs_lyr_src:
 525        if old_axis_mapping_srs is None and epsg_code_dest == 4326:
 526            old_axis_mapping_srs = drvr_name in DRIVERS_OLD_SRS_AXIS_MAPPING_4326
 527
 528        if epsg_code_dest:
 529            srs_lyr_dest = srs_ref_from_epsg_code(epsg_code_dest, old_axis_mapping=old_axis_mapping_srs)
 530        else:
 531            srs_lyr_dest = srs_lyr_src
 532
 533    if ds_gdal_dest.TestCapability(ODsCCreateLayer):
 534        # To avoid message "Warning"
 535        if ds_gdal_dest.TestCapability(ODsCDeleteLayer) and ds_gdal_dest.GetLayerByName(nom_layer):
 536            ds_gdal_dest.DeleteLayer(nom_layer)
 537
 538        layer_out, nom_layer = create_layer_on_ds_gdal(ds_gdal_dest, nom_layer, nom_geom_out, gtype, srs_lyr_dest,
 539                                                       **extra_opt_list)
 540    else:
 541        layer_out = ds_gdal_dest.GetLayer(0)
 542
 543    geom_field_out = None
 544    if layer_out.TestCapability(OLCAlterFieldDefn):
 545        layer_out_def = layer_out.GetLayerDefn()
 546        geom_field_out = get_geom_field_defn(layer_out_def)
 547        if geom_field_out:
 548            nom_geom_field_ref = act_geom_field.GetNameRef() if act_geom_field else None
 549            if not nom_geom_out:
 550                nom_geom_out = nom_geom_field_ref
 551                if not nom_geom_out:
 552                    nom_geom_out = "GEOMETRY"
 553            if not nom_geom_field_ref or nom_geom_field_ref.lower() != nom_geom_out.lower():
 554                geom_field_out.SetName(nom_geom_out)
 555
 556    if layer_out.TestCapability(OLCCreateField):
 557        for fd in fields_layer_gdal(layer_src):
 558            nom_fd = fd.GetNameRef().upper()
 559            if layer_out.FindFieldIndex(nom_fd, True) < 0 and (not sel_camps or nom_fd in sel_camps) and \
 560                    nom_fd not in (gn.upper() for gn in geoms_src_minus_affix):
 561                layer_out.CreateField(fd)
 562
 563        for gfd in geom_fields_layer_gdal(layer_src):
 564            nom_gfd = fix_affix_geom_name_layer_gdal(gfd.GetNameRef(), layer_src).upper()
 565            if (not sel_camps or nom_gfd in sel_camps) and \
 566                    (nom_gfd not in geoms_src or not exclude_cols_geoms) and \
 567                    nom_gfd != nom_geom_out:
 568                gfd_def_src = layer_src_def.GetGeomFieldDefn(layer_src_def.GetGeomFieldIndex(gfd.GetNameRef()))
 569                gfd_def_dest = GeomFieldDefn(nom_gfd, gfd_def_src.GetType())
 570                if srs_lyr_dest:
 571                    gfd_def_dest.SetSpatialRef(srs_lyr_dest)
 572                layer_out.CreateGeomField(gfd_def_dest)
 573
 574    # Si capa ALFA se fuerza a crear la capa aunque no tenga geometría
 575    if not nom_geom_out and not nom_geom_src:
 576        null_geoms = True  # Si no hay geom siempre se añaden
 577
 578    add_layer_features_to_layer(layer_src, ds_gdal_dest, layer_out, nom_geom_src, nom_geom_out,
 579                                srs_lyr_dest=srs_lyr_dest,
 580                                null_geoms=null_geoms,
 581                                tolerance_simplify=tolerance_simplify,
 582                                epsg_code_src_default=epsg_code_src_default,
 583                                old_axis_mapping_srs=old_axis_mapping_srs)
 584
 585    if drvr_name == DRVR_GPKG:
 586        create_spatial_index_layer_gpkg(ds_gdal_dest, nom_layer)
 587
 588    if drvr_name == DRVR_CSV:
 589        path_csv = ds_gdal_dest.GetName()
 590        if os.path.exists(path_csv):
 591            rename_wkt_geoms_csv(path_csv)
 592
 593    return ds_gdal_dest.GetLayerByName(nom_layer)
 594
 595
 596def rename_wkt_geoms_csv(path_csv):
 597    """
 598    Remove "_WKT" affix from WKT geometry fields in CSV file
 599
 600    Args:
 601        path_csv (str): Path to CSV file
 602
 603    Returns:
 604
 605    """
 606    with open(path_csv, 'r', encoding='utf-8') as r_file:
 607        content_csv = r_file.readlines()
 608        first_line = next(iter(content_csv), None)
 609    if first_line:
 610        new_first_line = first_line.replace(PREFFIX_GEOMS_LAYERS_GDAL_CSV, '')
 611        if new_first_line != first_line:
 612            content_csv[0] = new_first_line
 613            with open(path_csv, 'w', encoding='utf-8') as w_file:
 614                w_file.writelines(content_csv)
 615
 616
 617def add_layer_features_to_layer(layer_src, ds_gdal_dest, layer_dest, nom_geom_src=None, nom_geom_dest=None,
 618                                srs_lyr_dest=None, null_geoms=False, tolerance_simplify=None,
 619                                remove_prev_features=False, epsg_code_src_default=4326, old_axis_mapping_srs=None):
 620    """
 621    From a layer of a dataset, add the features to another layer of another dataset.
 622
 623    Args:
 624        layer_src (ogr.Layer): Layer origen
 625        ds_gdal_dest (gdal.Dataset): Dataset destino
 626        layer_dest (ogr.Layer): Layer destino
 627        nom_geom_src (str=None): Si el nombre no se corresponde con ninguna de las geometrias de la layer_src se utilizará
 628                            como ALIAS de la geometria 0 (la defecto) de la nueva layer resultante
 629        nom_geom_dest (str):
 630        srs_lyr_dest (ogr.SpatialReference=None): Spatial Reference System de la layer_dest
 631        null_geoms (bool=False): Por defecto no grabará las filas que la geometria principal (nom_geom) es nula
 632        tolerance_simplify (float=None): Tolerancia (distancia minima) en unidades del srs de la layer
 633             Mirar método Simplify() sobre osgeo.ogr.Geometry
 634        remove_prev_features (bool=False): Si es True, se eliminarán las features previas de la layer_dest
 635        epsg_code_src_default (int=4326): EPSG code del srs de la layer_src si no se encuentra
 636        old_axis_mapping_srs (bool=None): setea si quieres usar old mapping axes en (LONG, LAT) (GDAL >3.2)
 637
 638    Returns:
 639
 640    """
 641    geoms_out = geoms_layer_gdal(layer_dest)
 642    cols_out = cols_layer_gdal(layer_dest)
 643
 644    geom_transform = None
 645    srs_lyr_src = srs_for_layer(layer_src, nom_geom_src)
 646    if not srs_lyr_src:
 647        srs_lyr_src = srs_ref_from_epsg_code(epsg_code_src_default, old_axis_mapping_srs)
 648
 649    if srs_lyr_dest is None:
 650        srs_lyr_dest = srs_for_layer(layer_dest, nom_geom_dest)
 651
 652    if srs_lyr_dest and not srs_lyr_src.IsSame(srs_lyr_dest):
 653        geom_transform = osr.CoordinateTransformation(srs_lyr_src, srs_lyr_dest)
 654
 655    ds_trans = ds_gdal_dest.TestCapability(ODsCTransactions)
 656
 657    if remove_prev_features:
 658        if ds_trans:
 659            layer_dest.StartTransaction()
 660        for feat in layer_dest:
 661            layer_dest.DeleteFeature(feat.GetFID())
 662        if ds_trans:
 663            layer_dest.CommitTransaction()
 664
 665    if ds_trans:
 666        layer_dest.StartTransaction()
 667
 668    i = 0
 669    for feat_src, geom_src, nt_src in feats_layer_gdal(layer_src, nom_geom_src):
 670        cols_out_chk = [col.upper() for col in cols_out.union(geoms_out)]
 671        vals_camps = {nc: val for nc, val in nt_src._asdict().items()
 672                      if nc.upper() in cols_out_chk}
 673        if null_geoms or geom_src:
 674            if nom_geom_dest:
 675                vals_camps[nom_geom_dest.upper()] = geom_src
 676
 677            add_feature_to_layer_gdal(layer_dest,
 678                                      tolerance_simplify=tolerance_simplify,
 679                                      geom_trans=geom_transform,
 680                                      **vals_camps)
 681
 682            if i > 0 and (i % 1000) == 0:
 683                print_debug("{} registres tractats...".format(str(i)))
 684            i += 1
 685    if ds_trans:
 686        layer_dest.CommitTransaction()
 687
 688
 689def create_layer_on_ds_gdal(ds_gdal_dest, nom_layer, nom_geom=None, gtype=None, srs_lyr_out=None, **extra_opt_list):
 690    """
 691
 692    Args:
 693        ds_gdal_dest (ogr.DataSource):
 694        nom_layer (str):
 695        nom_geom (str=None):
 696        gtype (OGRwkbGeometryType=None): Indicar integer representativo del tipo de geometria de la layer (ogr.wkbPoint, ogr.wkbPolygon, ogr.LineString, ...)
 697        srs_lyr_out (osr.SpatialReference=None): codigo epsg del sistema de coordenadas de la geometria
 698        **extra_opt_list: Calves valores option list creacion layer gdal
 699
 700    Returns:
 701        layer_gdal (ogr.Layer), nom_layer (str)
 702    """
 703    drvr_name = ds_gdal_dest.GetDriver().GetName().upper()
 704
 705    opt_list = {k.upper(): v.upper() for k, v in extra_opt_list.items()}
 706
 707    opt_geom_name = "GEOMETRY_NAME"
 708    if nom_geom:
 709        opt_list[opt_geom_name] = "{}={}".format(opt_geom_name, nom_geom.upper())
 710    else:
 711        if opt_geom_name in opt_list:
 712            opt_list.pop(opt_geom_name)
 713        gtype = ogr.wkbNone
 714
 715    opt_list["IDENTIFIER"] = opt_list.get("IDENTIFIER", "IDENTIFIER={}".format(nom_layer))
 716
 717    opt_list = set_create_option_list_for_driver_gdal(drvr_name, **{k.upper(): v.upper() for k, v in opt_list.items()})
 718
 719    if ds_gdal_dest.TestCapability(ODsCDeleteLayer) and ds_gdal_dest.GetLayerByName(nom_layer):
 720        print_warning("!ATENCION! - Se sobreescribirá la layer '{}' sobre el datasource GDAL '{}'".format(
 721            nom_layer, ds_gdal_dest.GetDescription()))
 722        ds_gdal_dest.DeleteLayer(nom_layer)
 723
 724    if drvr_name.upper() == DRVR_KML:
 725        nom_layer = nom_layer.replace("-", "__")
 726
 727    layer_out = ds_gdal_dest.CreateLayer(nom_layer, srs_lyr_out, geom_type=gtype, options=list(opt_list.values()))
 728
 729    return layer_out, layer_out.GetName() if layer_out else None
 730
 731
 732def create_spatial_index_layer_gpkg(ds_gpkg, nom_layer):
 733    """
 734    Crea spatial index sobre una layer (layer_gpkg) de un datosource gpkg (ds_gpkg)
 735
 736    Args:
 737        ds_gpkg (osgeo.ogr.DataSource):
 738        layer_gpkg (ogr.Layer):
 739
 740    Returns:
 741        bool
 742    """
 743    layer_gpkg = ds_gpkg.GetLayerByName(nom_layer)
 744    if layer_gpkg and (geom_name := layer_gpkg.GetGeometryColumn()):
 745        tab_name = layer_gpkg.GetName()
 746        res = ds_gpkg.ExecuteSQL(f"SELECT HasSpatialIndex('{tab_name}', '{geom_name}');")
 747        if not bool(res and res.GetNextFeature().GetField(0) == 1):
 748            # Se crea spatial_index ya que GDAL NO lo hace
 749            ds_gpkg.StartTransaction()
 750            ds_gpkg.ExecuteSQL(f"SELECT CreateSpatialIndex('{tab_name}', '{geom_name}');")
 751            ds_gpkg.CommitTransaction()
 752
 753        return True
 754    else:
 755        return False
 756
 757
 758def add_feature_to_layer_gdal(layer_gdal, tolerance_simplify=None, geom_trans=None, commit=False, **valors_camps):
 759    """
 760
 761    Args:
 762        layer_gdal (ogr.Layer):
 763        tolerance_simplify (float=None): Tolerancia (distancia minima) en unidades del srs de la layer.
 764                                        Mirar método Simplify() sobre osgeo.ogr.Geometry
 765        geom_trans (osr.CoordinateTransformation=None): transformacion para convertir las geometrias a otro SRS
 766        commit (bool=False): Per defecte no farà commit de la transaccio
 767        **valors_camps: pares nombre_campo=valor de la feature a crear
 768
 769    Returns:
 770        feat (ogr.Feature)
 771    """
 772    dd_feat = ogr.Feature(layer_gdal.GetLayerDefn())
 773
 774    for camp, val in valors_camps.items():
 775        idx_geom = dd_feat.GetGeomFieldIndex(camp)
 776        es_geom = idx_geom >= 0
 777        if es_geom:
 778            if val:
 779                new_val = None
 780                if tolerance_simplify:
 781                    new_val = val.Simplify(tolerance_simplify)
 782                if val and geom_trans:
 783                    new_val = val.Transform(geom_trans)
 784                val = new_val or val
 785
 786            dd_feat.SetGeomField(idx_geom, val)
 787
 788        idx_fld = dd_feat.GetFieldIndex(camp)
 789        if idx_fld >= 0:
 790            if val:
 791                if es_geom or hasattr(val, "ExportToIsoWkt"):
 792                    val = val.ExportToIsoWkt()
 793                dd_feat.SetField(camp, val)
 794            else:
 795                dd_feat.SetFieldNull(camp)
 796
 797        if idx_fld < 0 and idx_geom < 0:
 798            if isinstance(val, Geometry):
 799                if geom_trans:
 800                    val.Transform(geom_trans)
 801                dd_feat.SetGeometry(val)
 802            elif val:
 803                print_warning("!ATENCION! - La :layer_gdal no contiene el campo '{}'".format(camp))
 804
 805    commit_trans = commit and layer_gdal.TestCapability(OLCTransactions)
 806    if commit_trans:
 807        layer_gdal.StartTransaction()
 808
 809    new_feat = layer_gdal.CreateFeature(dd_feat)
 810
 811    if commit_trans:
 812        layer_gdal.CommitTransaction()
 813
 814    return new_feat
 815
 816
 817def drivers_ogr_gdal_disponibles():
 818    """
 819    Retorna lista de drivers disponibles a través de la librería osgeo-gdal-ogr
 820
 821    Returns:
 822        dict
 823    """
 824    cnt = ogr.GetDriverCount()
 825    driver_list = []
 826    drivers = OrderedDict()
 827
 828    for i in range(cnt):
 829        driver = ogr.GetDriver(i)
 830        driver_name = driver.GetName()
 831        driver_list.append(driver_name)
 832
 833    for driver_name in driver_list:
 834        # Is File GeoDatabase available?
 835        drv = ogr.GetDriverByName(driver_name)
 836        if drv is None:
 837            print_warning("{} !!ATENTION - driver NOT available!!".format(driver_name))
 838        else:
 839            drivers[driver_name] = drv
 840            print_debug(driver_name)
 841
 842    return drivers
 843
 844
 845def drivers_ogr_gdal_vector_file():
 846    """
 847    Devuelve diccionario con los driver gdal para fichero vectorial
 848
 849    Returns:
 850        dict
 851    """
 852    result = {}
 853    for nd, d in drivers_ogr_gdal_disponibles().items():
 854        try:
 855            # Try GetMetadata_Dict first (older GDAL versions)
 856            if hasattr(d, "GetMetadata_Dict"):
 857                metadata = d.GetMetadata_Dict()
 858            else:
 859                # Fall back to GetMetadata for GDAL 3.10+
 860                metadata = d.GetMetadata()
 861
 862            if metadata and metadata.get('DMD_EXTENSIONS'):
 863                result[nd] = d
 864        except Exception:
 865            # Skip drivers that have issues
 866            continue
 867
 868    return result
 869
 870
 871def format_nom_column(nom_col, case_format=None):
 872    """
 873    Retorna el nombre de columna formateado para usar como campo de un namedtuple (quita espacios y los reemplaza por "_").
 874
 875    Args:
 876        nom_col (str):
 877        case_format (str|None): formato de la columna. Si 'upper' or 'lower'. Por defecto tal como viene indicado
 878
 879    Returns:
 880        str
 881    """
 882    nom_col_res = nom_col.replace(" ", "_")
 883    if case_format:
 884        if (case_format := case_format.lower()) == 'upper':
 885            nom_col_res = nom_col_res.upper()
 886        elif case_format == 'lower':
 887            nom_col_res = nom_col_res.lower()
 888
 889    return nom_col_res
 890
 891
 892def namedtuple_layer_gdal(layer_gdal, extract_affix_geom_fld=None, case_fields_nt=None):
 893    """
 894    Devuelve namedTuple con los campos del layer pasado por parametro
 895
 896    Args:
 897        layer_gdal:
 898        extract_affix_geom_fld (str=None):
 899        case_fields_nt (str|None): formato de las columnas de NamedTuple retornado. Si 'upper' or 'lower'. Por defecto tal como vienen de la FeatureLayer
 900
 901    Returns:
 902        namedtuple: con nombre "gdalFeatDef_{NOM_LAYER}" y con los campos de la layer
 903    """
 904    camps_layer = set()
 905    for fld in fields_layer_gdal(layer_gdal):
 906        camps_layer.add(format_nom_column(fld.GetNameRef(), case_format=case_fields_nt))
 907
 908    for nom_geom in geoms_layer_gdal(layer_gdal, extract_affix=extract_affix_geom_fld):
 909        camps_layer.add(format_nom_column(nom_geom, case_format=case_fields_nt))
 910
 911    nom_layer = layer_gdal.GetName().upper().split(".")[0].replace("-", "_")
 912    return namedtuple(f"gdalFeatDef_{nom_layer}", camps_layer)
 913
 914
 915def feats_layer_ds_gdal(ds_gdal, nom_layer=None, filter_sql=None, case_fields_nt=None):
 916    """
 917    Itera las features (registros de una layer de gdal) y los devuelve como un namdetuple
 918
 919    Args:
 920        ds_gdal: datasource gdal
 921        nom_layer (str=None): Si no viene informado cogerá la primera layer que encuentre en el datasource
 922        filter_sql (str=None): Si viene informado se aplicará como filtro sql a la layer seleccionada.
 923                            Utiliza OGR SQL (vease https://www.gdal.org/ogr_sql.html)
 924        case_fields_nt (str|None): formato de las columnas de NamedTuple retornado. Si 'upper' or 'lower'. Por defecto tal como vienen de la FeatureLayer
 925
 926    Returns:
 927        ogr.Feature, ogr.Geometry, namedtuple_layer_gdal
 928    """
 929    if not nom_layer:
 930        layer_gdal = ds_gdal.GetLayer()
 931    else:
 932        layer_gdal = ds_gdal.GetLayerByName(nom_layer)
 933
 934    if layer_gdal:
 935        for feat, geom, vals in feats_layer_gdal(layer_gdal, filter_sql=filter_sql, case_fields_nt=case_fields_nt):
 936            yield feat, geom, vals
 937    else:
 938        raise ValueError(f"'{nom_layer}' is not a valid layer over DataSource GDAL '{ds_gdal.GetName()}'")
 939
 940
 941def feats_layer_gdal(layer_gdal, nom_geom=None, filter_sql=None, extract_affix_geom_fld=PREFFIX_GEOMS_LAYERS_GDAL,
 942                     case_fields_nt=None):
 943    """
 944    Itera las features (registros de una layer de gdal) y los devuelve como un namedtuple
 945
 946    Args:
 947        layer_gdal (ogr.Layer):
 948        nom_geom (str=None): Por defecto la geometria activa o principal
 949        filter_sql (str=None): Si viene informado se aplicará como filtro sql a la layer seleccionada.
 950                            Utiliza OGR SQL (vease https://www.gdal.org/ogr_sql.html)
 951        extract_affix_geom_fld (str=SUFFIX_GEOMS_LAYERS_GDAL): Por defecto quita el sufijo 'geom_' a los campos geom
 952        case_fields_nt (str=None): si 'upper' los campos del NamedTuple retornado estarán en UPPERCASE y si 'lower' al contrario.
 953                        Por defecto tal como vienen de la FeatureLayer
 954
 955    Returns:
 956        ogr.Feature, ogr.Geometry, namedtuple_layer_gdal
 957    """
 958    layer_gdal.ResetReading()
 959    ntup_layer = namedtuple_layer_gdal(layer_gdal, extract_affix_geom_fld, case_fields_nt=case_fields_nt)
 960    n_geoms_layer = {nom_geom_feat: fix_affix_geom_name_layer_gdal(nom_geom_feat, layer_gdal, extract_affix_geom_fld)
 961                     for nom_geom_feat in geoms_layer_gdal(layer_gdal)}
 962
 963    if filter_sql:
 964        layer_gdal.SetAttributeFilter(filter_sql)
 965
 966    def vals_feature_gdal(feat_gdal):
 967        vals = {}
 968        feat_items = feat_gdal.items()
 969        # In GDAL 3.x, items() returns a dict
 970        if isinstance(feat_items, dict):
 971            for camp, val in feat_items.items():
 972                vals[format_nom_column(camp, case_format=case_fields_nt)] = val
 973        else:
 974            # Older GDAL versions may return different structure
 975            for camp, val in feat_items:
 976                vals[format_nom_column(camp, case_format=case_fields_nt)] = val
 977
 978        for camp_geom, camp_geom_val in n_geoms_layer.items():
 979            idx_geom = feat_gdal.GetGeomFieldIndex(camp_geom)
 980            if idx_geom >= 0:
 981                vals[format_nom_column(camp_geom_val, case_format=case_fields_nt)] = feat_gdal.GetGeomFieldRef(idx_geom)
 982
 983        return vals
 984
 985    if layer_gdal:
 986        for f_tab in layer_gdal:
 987            idx_geom = f_tab.GetGeomFieldIndex(nom_geom) if nom_geom else -1
 988            yield f_tab, \
 989                f_tab.GetGeomFieldRef(idx_geom) if idx_geom >= 0 else f_tab.geometry(), \
 990                ntup_layer(**vals_feature_gdal(f_tab))
 991
 992        layer_gdal.ResetReading()
 993
 994
 995def distinct_vals_camp_layer_gdal(layer_gdal, nom_camp, filter_sql=None):
 996    """
 997    Devuelve set con distintos valores para el campo indicado del layer GDAL indicada
 998
 999    Args:
1000        layer_gdal (ogr.Layer):
1001        nom_camp (str):
1002        filter_sql (str=None):
1003
1004    Returns:
1005        set
1006    """
1007    if nom_camp.upper() not in cols_layer_gdal(layer_gdal):
1008        raise Exception("Argumento :NOM_CAMP = '{}' erróneo ya que "
1009                        "LAYER_GDAL no contiene el campo indicado".format(nom_camp))
1010
1011    return {getattr(nt_feat, nom_camp.upper())
1012            for feat, geom, nt_feat in feats_layer_gdal(layer_gdal, filter_sql=filter_sql)}
1013
1014
1015def fields_layer_gdal(layer_gdal):
1016    """
1017    Itera sobre los FieldDefn de una layer gdal
1018
1019    Args:
1020        layer_gdal:
1021
1022    Yields:
1023        osgeo.ogr.FieldDefn
1024    """
1025    layer_def = layer_gdal.GetLayerDefn()
1026    for i in range(0, layer_def.GetFieldCount()):
1027        yield layer_def.GetFieldDefn(i)
1028
1029    layer_def = None
1030
1031
1032def geom_fields_layer_gdal(layer_gdal):
1033    """
1034    Itera sobre los GeomFieldDefn de una layer gdal
1035
1036    Args:
1037        layer_gdal:
1038
1039    Yields:
1040        osgeo.ogr.GeomFieldDefn
1041    """
1042    layer_def = layer_gdal.GetLayerDefn()
1043    for i in range(0, layer_def.GetGeomFieldCount()):
1044        yield layer_def.GetGeomFieldDefn(i)
1045
1046    layer_def = None
1047
1048
1049def nom_layers_datasource_gdal(ds_gdal):
1050    """
1051    Retorna lista con los nombres de las layers de un datasource gdal.
1052    CACHEA esos nombres por DS Name para optimizar llamadas repetidas.
1053
1054    Args:
1055        ds_gdal (ogr.Datasource):
1056
1057    Returns:
1058        set
1059    """
1060    if not ds_gdal:
1061        return set()
1062
1063    ds_name = ds_gdal.GetName() if hasattr(ds_gdal, "GetName") else None
1064    layer_count = ds_gdal.GetLayerCount() if hasattr(ds_gdal, "GetLayerCount") else None
1065
1066    if ds_name and ds_name in _CACHE_NOM_LAYERS_DS_GDAL:
1067        cached = _CACHE_NOM_LAYERS_DS_GDAL[ds_name]
1068        if layer_count is None or cached.get("count") == layer_count:
1069            # Devolvemos copia para evitar mutaciones externas de la caché
1070            return set(cached.get("layers", set()))
1071
1072    layers = {l.GetName() for l in ds_gdal}
1073
1074    # Solo cacheamos cuando existe identificador estable de datasource
1075    if ds_name:
1076        _CACHE_NOM_LAYERS_DS_GDAL[ds_name] = {
1077            "count": layer_count,
1078            "layers": set(layers)
1079        }
1080
1081    return layers
1082
1083
1084def cols_layer_gdal(layer_gdal):
1085    """
1086    Retorna lista con las columnas de una layer gdal
1087
1088    Args:
1089        layer_gdal:
1090
1091    Returns:
1092        set
1093    """
1094    camps = set()
1095    for fd in fields_layer_gdal(layer_gdal):
1096        # camps.add(fd.GetName().upper())
1097        camps.add(fd.GetNameRef())
1098
1099    return camps
1100
1101
1102def geoms_layer_gdal(layer_gdal, extract_affix=None):
1103    """
1104    Retorna lista con las columnas geométricas de una layer gdal
1105
1106    Args:
1107        layer_gdal:
1108        extract_affix (str=None=): if informed extract the affix passed
1109
1110    Returns:
1111        set
1112    """
1113    camps_geom = set()
1114    for gdf in geom_fields_layer_gdal(layer_gdal):
1115        # camps_geom.add(gdf.GetName().upper())
1116        name_ref = gdf.GetNameRef()
1117        if extract_affix:
1118            name_ref = fix_affix_geom_name_layer_gdal(name_ref, layer_gdal, extract_affix)
1119        camps_geom.add(name_ref)
1120
1121    return camps_geom
1122
1123
1124def fix_affix_geom_name_layer_gdal(geom_name, layer_gdal, affix=PREFFIX_GEOMS_LAYERS_GDAL):
1125    """
1126    Extract affix (default=PREFFIX_GEOMS_LAYERS_GDAL) from name geoms if correspond with other column
1127    Fix GDAL > 3.4 where geoms on GeoCSV arrive with preffix 'geom_'
1128
1129    Args:
1130          geom_name (str):
1131          layer_gdal:
1132          affix (str=PREFFIX_GEOMS_LAYERS_GDAL):
1133
1134    Returns:
1135        str
1136    """
1137    cols_layer = cols_layer_gdal(layer_gdal)
1138    geom_name_layer = geom_name
1139    if geom_name.lower().startswith(affix):
1140        geom_name_aux = re.sub(affix, '', geom_name, flags=re.IGNORECASE)
1141        if geom_name_aux.lower() in (col.lower() for col in cols_layer):
1142            geom_name_layer = geom_name_aux
1143
1144    return geom_name_layer
1145
1146
1147def add_layer_gdal_to_ds_gdal(ds_gdal, layer_gdal, nom_layer=None, lite=False, srs_epsg_code=None, multi_geom=False,
1148                              nom_geom=None, null_geoms=False, **extra_opt_list):
1149    """
1150    Añade una layer_gdal a un datasource_gdal. Si es una layer con multigeometrias las separa en una layer por geometria
1151
1152    Args:
1153        ds_gdal (osgeo.ogr.Datasource):
1154        layer_gdal (osgeo.ogr.Layer):
1155        nom_layer (str=None):
1156        lite (bool=False):
1157        srs_epsg_code (int=None): codigo EPSG para el sistema de coordenadas con el que se quieren convertir
1158                                  las geometrias
1159        nom_geom (str=None): nombre de geometria de layer origen (layer_gdal) que se copiará, si no todas
1160        multi_geom (bool=False): Si el DS_GDAL destino permite multigeometria
1161        null_geoms (bool=False): Indica si se admitirán registros con NULL geoms. Por defecto NO
1162        **extra_opt_list (str): Lista claves-valores de opciones para createLayer del driver de ds_gdal indicado
1163
1164    Returns:
1165        new_layers_ds_gdal (list)
1166    """
1167    new_layers_ds_gdal = []
1168
1169    if not nom_layer:
1170        nom_layer = layer_gdal.GetName()
1171
1172    geoms_layer = geoms_layer_gdal(layer_gdal)
1173    if nom_geom:
1174        if nom_geom.upper() not in (g.upper() for g in geoms_layer):
1175            Exception("!ERROR! - Nombre de geometria '{}' no existe en la layer GDAL origen")
1176        else:
1177            geoms_layer = (nom_geom,)
1178
1179    if geoms_layer:
1180        tol = None
1181        if lite:
1182            tol = SIMPLIFY_TOLERANCE
1183
1184        nom_layer_base = nom_layer.split("-")[0]
1185        if not multi_geom:
1186            for geom_name in geoms_layer:
1187                # Fix GDAL > 3.4 where geoms on GeoCSV arrive with preffix 'geom_'
1188                geom_name_layer = fix_affix_geom_name_layer_gdal(geom_name, layer_gdal)
1189
1190                nom_layer = "{}-{}".format(nom_layer_base, geom_name_layer).lower()
1191                extra_opt_list["GEOMETRY_NAME"] = "GEOMETRY_NAME={}".format(geom_name_layer)
1192                lyr = create_layer_from_layer_gdal_on_ds_gdal(ds_gdal, layer_gdal, nom_layer, geom_name,
1193                                                              tolerance_simplify=tol, null_geoms=null_geoms,
1194                                                              epsg_code_dest=srs_epsg_code, **extra_opt_list)
1195                new_layers_ds_gdal.append(lyr)
1196        else:
1197            lyr = create_layer_from_layer_gdal_on_ds_gdal(ds_gdal, layer_gdal, nom_layer, exclude_cols_geoms=False,
1198                                                          tolerance_simplify=tol, null_geoms=True,
1199                                                          epsg_code_dest=srs_epsg_code, **extra_opt_list)
1200            new_layers_ds_gdal.append(lyr)
1201    else:
1202        if ds_gdal.GetDriver().GetName() == DRVR_PG:
1203            extra_opt_list.update(dict(
1204                precision=extra_opt_list.get('precision', 'PRECISION=NO')
1205            ))
1206
1207        lyr = copy_layer_gdal_to_ds_gdal(layer_gdal, ds_gdal, nom_layer.lower(), **extra_opt_list)
1208        new_layers_ds_gdal.append(lyr)
1209
1210    return new_layers_ds_gdal
1211
1212
1213def copy_layers_gpkg(ds_gpkg, driver, dir_base, lite=False, srs_epsg_code=None, zipped=True):
1214    """
1215
1216    Args:
1217        ds_gpkg (osgeo.ogr.Datasource):
1218        driver (str):
1219        dir_base (str=None):
1220        lite (bool=False):
1221        srs_epsg_code (int=None): codigo EPSG para el sistema de coordenadas con el que se quieren convertir
1222                                  las geometrias
1223        zipped (bool=False):
1224
1225    Returns:
1226        num_layers (int)
1227    """
1228    num_layers = 0
1229    subdir_drvr = os.path.normpath(os.path.join(dir_base, driver.upper()))
1230    utils.create_dir(subdir_drvr)
1231
1232    for layer_gpkg in (ds_gpkg.GetLayer(id_lyr) for id_lyr in range(ds_gpkg.GetLayerCount() - 1)):
1233        if driver == DRVR_GPKG:
1234            nom_ds, ext = utils.split_ext_file(os.path.basename(ds_gpkg.GetName()))
1235        else:
1236            nom_ds = f"{layer_gpkg.GetName()}".lower()
1237
1238        ds_gdal, existia = datasource_gdal_vector_file(driver, nom_ds, subdir_drvr)
1239
1240        add_layer_gdal_to_ds_gdal(ds_gdal, layer_gpkg, lite=lite, srs_epsg_code=srs_epsg_code)
1241
1242        num_layers += 1
1243
1244    return num_layers
1245
1246
1247def set_csvt_for_layer_csv(path_csv, **tipus_camps):
1248    """
1249    Crea/Modifica el CSVT asociado con los tipos indicados para cada columna
1250
1251    Args:
1252        path_csv (str):
1253        **tipus_camps: clave=valor con el nombre del campo y el tipo de campo asociado (p.e. String(25), WKT, Integer,...)
1254
1255    Returns:
1256        path_csvt (str)
1257    """
1258    lyr_csv, nom_layer, ds_lyr = layer_gdal_from_file(path_csv, DRVR_CSV)
1259    if not lyr_csv:
1260        print_warning("!ATENCIO! - No s'ha pogut obrir la layer CSV '{}'".format(path_csv))
1261        return
1262
1263    path_lyr_csvt = os.path.join(os.path.dirname(path_csv), "{}.csvt".format(nom_layer))
1264    tips_lyr = {}
1265    for fld in fields_layer_gdal(lyr_csv):
1266        tip_fld = fld.GetFieldTypeName(fld.GetType())
1267        sufix = ""
1268        w = fld.GetWidth()
1269        if w:
1270            sufix = "{}".format(w)
1271        p = fld.GetPrecision()
1272        if p:
1273            sufix += ".{}".format(p)
1274        if sufix:
1275            tip_fld += "({})".format(sufix)
1276
1277        tips_lyr[fld.name.upper()] = tip_fld
1278
1279    for nom_camp, tip_camp in tipus_camps.items():
1280        nom_camp = nom_camp.upper()
1281        if nom_camp not in tips_lyr:
1282            print_warning("!ATENCIO! - Camp '{}' no existeix sobre la layer CSV '{}'".format(nom_camp, path_csv))
1283            continue
1284
1285        tips_lyr[nom_camp] = tip_camp
1286
1287    with open(path_lyr_csvt, mode="w", encoding="utf8") as f_csvt:
1288        f_csvt.write(",".join(tips_lyr.values()))
1289
1290
1291def zip_and_clean_ds_csv(path_csv):
1292    """
1293    # Zipea datasource osgeo CSV y como GDAL no crea los tipos WKT para las geometrias en el CSVT se fuerza
1294
1295    Args:
1296        path_csv (str): path Datasource osgeo CSV
1297
1298    Returns:
1299        zip_path (str)
1300    """
1301    ds_gdal_csv, dummy = datasource_gdal_vector_file(DRVR_CSV, os.path.basename(path_csv).split(".")[0],
1302                                                     os.path.dirname(path_csv))
1303    layer_csv = ds_gdal_csv.GetLayer(0)
1304    tips_geoms = {gn: "WKT" for gn in
1305                  (*geoms_layer_gdal(layer_csv, extract_affix=PREFFIX_GEOMS_LAYERS_GDAL),
1306                   *geoms_layer_gdal(layer_csv, extract_affix=PREFFIX_GEOMS_LAYERS_GDAL_CSV))}
1307    path_csv = ds_gdal_csv.GetDescription()
1308    ds_gdal_csv = None  # Cerramos el datasource para que se guarden los cambios
1309    if path_csv and os.path.exists(path_csv):
1310        if tips_geoms:
1311            set_csvt_for_layer_csv(path_csv, **tips_geoms)
1312        rename_wkt_geoms_csv(path_csv)
1313    dir_base_csv = os.path.dirname(path_csv)
1314    nom_layer = os.path.basename(path_csv).split(".")[0]
1315    l_files_csv = [os.path.join(dir_base_csv, nf) for nf in os.listdir(dir_base_csv)
1316                   if nf.lower().startswith(nom_layer.lower()) and
1317                   any(nf.lower().endswith(ext) for ext in ('csv', 'csvt'))]
1318    zip_path = utils.zip_files(os.path.join(dir_base_csv, "{}.zip".format(nom_layer)), l_files_csv)
1319    if zip_path:
1320        for fl_csv in l_files_csv:
1321            os.remove(fl_csv)
1322
1323    return zip_path
1324
1325
1326def convert_angle(pt_xy, deg_ang, orig_srs, dest_srs):
1327    """
1328
1329    Args:
1330        pt_xy (tuple):
1331        deg_ang (float):
1332        orig_srs (osr.SpatialReference):
1333        dest_srs (osr.SpatialReference):
1334
1335    Returns:
1336
1337    """
1338    trans = osr.CoordinateTransformation(orig_srs, dest_srs)
1339    x, y = pt_xy
1340    dx = math.sin(math.radians(deg_ang)) * 0.00000001
1341    dy = math.cos(math.radians(deg_ang)) * 0.00000001
1342    pt1 = ogr.CreateGeometryFromWkt("POINT ({} {})".format(x, y))
1343    pt2 = ogr.CreateGeometryFromWkt("POINT ({} {})".format(x + dx, y + dy))
1344    pt1.Transform(trans)
1345    pt2.Transform(trans)
1346    x1, y1, z1 = pt1.GetPoint()
1347    x2, y2, z2 = pt2.GetPoint()
1348
1349    return math.degrees(math.atan2(y2 - y1, x2 - x1))
1350
1351
1352def transform_ogr_geom(a_ogr_geom, from_espg_code, to_epsg_code):
1353    """
1354    Transforma una geometria OGR según los EPSG indicados
1355
1356    Args:
1357        a_ogr_geom (ogr.geometry): una geometria del tipo OGR
1358        from_espg_code (int): codigo numérico del EPSG actual para la geometria
1359        to_epsg_code (int): codigo numérico del EPSG al que se quiere transformar
1360
1361    Returns:
1362        ogr.geometry
1363    """
1364    source = osr.SpatialReference()
1365    source.ImportFromEPSG(from_espg_code)
1366
1367    target = osr.SpatialReference()
1368    target.ImportFromEPSG(to_epsg_code)
1369
1370    a_transform = osr.CoordinateTransformation(source, target)
1371    a_ogr_geom.Transform(a_transform)
1372
1373    return a_ogr_geom
1374
1375
1376def ds_postgis(pg_conn_str, schemas=None, **open_options):
1377    """
1378    Retorna datasource GDAL para ddbb postgis
1379
1380    Args:
1381        pg_conn_str: str "PG:user=... password=..."
1382        schemas:
1383        **open_options: Open options del driver PostgreSQL de GDAL
1384                        (vease https://gdal.org/en/stable/drivers/vector/pg.html#dataset-open-options)
1385                        p.e.: LIST_ALL_TABLES='YES', SKIP_VIEWS='YES', SCHEMAS='public,myschema',
1386                              TABLES='tabla1,tabla2', OGR_PG_RETRIEVE_FID='YES'
1387
1388    Returns:
1389        osgeo.ogr.DataSource
1390    """
1391    if schemas:
1392        if isinstance(schemas, str):
1393            schemas = [*schemas.split(",")]
1394
1395        active_schema = schemas[0]
1396        str_schemas = ",".join(schemas)
1397        pg_conn_str = f"{pg_conn_str} active_schema='{active_schema}' schemas='{str_schemas}'"
1398
1399    if open_options:
1400        oo_list = [f"{k.upper()}={v}" for k, v in open_options.items()]
1401        return gdal.OpenEx(pg_conn_str, gdal.OF_VECTOR | gdal.OF_UPDATE, open_options=oo_list)
1402
1403    drvr, exts = driver_gdal(DRVR_PG)
1404
1405    return drvr.Open(pg_conn_str, 1)
1406
1407
1408def pg_string_connection(dbname=None, host=None, port=None, service=None, user=None, password=None, schemas=None):
1409    """
1410    Retorna una cadena de conexion a una base de datos postgis
1411
1412    Args:
1413        dbname:
1414        host:
1415        port:
1416        service:
1417        user:
1418        password:
1419        schemas (str | list[str]=None): Esquema de la base de datos
1420
1421    Returns:
1422        str
1423    """
1424    pg_conn = "PG:"
1425    if dbname:
1426        pg_conn += f"dbname='{dbname}'"
1427    if host:
1428        pg_conn += f"host='{host}'"
1429    if port:
1430        pg_conn += f"port='{port}'"
1431    if service:
1432        pg_conn += f"service='{service}'"
1433    if user:
1434        pg_conn += f"user='{user}'"
1435    if password:
1436        pg_conn += f"password='{password}'"
1437
1438    if schemas:
1439        if isinstance(schemas, str):
1440            schemas = [*schemas.split(",")]
1441
1442        active_schema = schemas[0]
1443        str_schemas = ",".join(schemas)
1444        pg_conn_str = f"{pg_conn} active_schema='{active_schema}' schemas='{str_schemas}'"
1445
1446    return pg_conn
1447
1448
1449def get_ds_postgis(dbname=None, host=None, port=None, service=None, user=None, password=None, schemas=None,
1450                   pg_conn_str=None, logger: Logger = None, **open_options) -> ogr.DataSource:
1451    """
1452    Retorna cached datasource postgis GDAL
1453
1454    Args:
1455        dbname (str=None): Nombre de la base de datos
1456        host (str=None): Host de la base de datos
1457        port (str=None): Puerto de la base de datos
1458        service (str=None): nombre service definido en pg_service.conf (env var PGSERVICE)
1459        user (str=None): Usuario de la base de datos
1460        password (str=None): Contraseña de la base de datos
1461        schemas (str | list[str]=None): Esquema de la base de datos
1462        pg_conn_str (str=None): cadena conexion tipo "PG:user=... password=... service=...". Si viene informada no se hace caso a ningún otro argumento
1463        logger (logging.Logger=None): (opcional) logger para LOG.
1464        **open_options (str): Opciones para abrir el Datasource de Postgis. Vease https://gdal.org/en/stable/drivers/vector/pg.html#dataset-open-options
1465    Returns:
1466        datasource_gdal (osgeo.gdal.DataSource)
1467    """
1468    ds_gdal: ogr.DataSource | None = None
1469
1470    if not pg_conn_str:
1471        pg_conn_str = pg_string_connection(dbname=dbname, host=host, port=port, service=service, user=user,
1472                                           password=password)
1473    elif not pg_conn_str.startswith('PG:'):
1474        pg_conn_str = f"PG:{pg_conn_str}"
1475
1476    # Retorna la pg_conn_str con el valor despues de 'password=' substituido por '****'
1477    pg_id_ds = re.sub(r"(?i)(password\s*=\s*')([^']*)(')", r"\1****\3", pg_conn_str)
1478    
1479    if ds_gdal := _CACHE_DS_POSTGIS.get(pg_id_ds):
1480        # Antes de retornar verificar que la conexión sigue siendo correcta
1481        try:
1482            _ = ds_gdal.GetLayerCount()
1483            lyr_test = ds_gdal.ExecuteSQL("SELECT 1")
1484            if lyr_test is not None:
1485                ds_gdal.ReleaseResultSet(lyr_test)
1486            return ds_gdal
1487        except Exception:
1488            _CACHE_DS_POSTGIS.pop(pg_id_ds, None)
1489
1490    if not logger:
1491        logger = utils_logging.get_base_logger()
1492
1493    error_exc = None
1494    try:
1495        ds_gdal: ogr.DataSource = ds_postgis(
1496            pg_conn_str=pg_conn_str,
1497            schemas=schemas,
1498            **open_options)
1499    except Exception as exc:
1500        error_exc = f"Error: {exc}"
1501    finally:
1502        if ds_gdal:
1503            logger.info(f"Conectado a la base de datos POSTGIS del repositorio GIS: {pg_id_ds}")
1504        else:
1505            msg_error = f"No se ha podido conectar a la base de datos POSTGIS [{pg_id_ds}] del repositorio GIS. " \
1506                        f"{error_exc if error_exc else ''} "
1507            logger.error(msg_error)
1508            raise ConnectionError(msg_error)
1509
1510    _CACHE_DS_POSTGIS[pg_id_ds] = ds_gdal
1511
1512    return ds_gdal
1513
1514
1515def reset_sequence_layer_postgis(ds_postgis, nom_layer):
1516    """
1517    Resetea la secuencia de la layer PostGIS indicada
1518
1519    Args:
1520        ds_postgis (osgeo.ogr.DataSource): datasource de la ddbb PostGIS:
1521        nom_layer (str): nombre de la layer PostGIS
1522
1523    Returns:
1524        bool: True si se ha reseteado correctamente
1525    """
1526    ok = False
1527    layer_dest = ds_postgis.GetLayerByName(nom_layer)
1528    fid_column = layer_dest.GetFIDColumn()
1529    get_seq_sql = f"SELECT pg_get_serial_sequence('{nom_layer}', '{fid_column}') FROM information_schema.columns " \
1530                  f"WHERE table_name = '{nom_layer}' AND column_default LIKE 'nextval%'"
1531
1532    if res_get_seq := ds_postgis.ExecuteSQL(get_seq_sql):
1533        seq_name = res_get_seq.GetNextFeature().GetField(0)
1534        ds_postgis.StartTransaction()
1535        reset_seq = f"SELECT setval('{seq_name}', 1, false);"
1536        ds_postgis.ExecuteSQL(reset_seq)
1537        ds_postgis.CommitTransaction()
1538        ok = True
1539
1540    return ok
1541
1542
1543if __name__ == '__main__':
1544    import fire
1545
1546    fire.Fire()