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