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