Etiqueta: python

  • 4.1 Rssi vs Distancia. Linealiza UN intervalo

    Referencia: Chapra 17.1 p 466. Burden 8.1 p498, Mínimos cuadrados en Métodos numéricos

    La linealización de curvas con Método de mínimos cuadrados se realiza usando las funciones de Numpy: np.linalg.lstsq()

    Para el análisis de una baliza, 'gtwFIEC', se obtienen los datos desde el archivo 'resumen_RssiUbica01' obtenido en:

    Integrar las tablas de Rssi y coordenadas de los puntos

    entregando como resultado un archivo con las ecuaciones obtenidas: 'resumen_ecuacionSimple05.json'

    La selección de la baliza  se realiza con un diccionario indicando la acción de 'analizar' como verdadero o falso (1,0), entre otros parámetros.

     'gtwFIEC':{'analizar'  : 1,
                'atipico_std' : 1,
                'grp' : ['FIEC','FCNM'],
                'tip' : ['punto'],
                'LOS' : [1,0]}
    

    Los valores atípicos se los discrimina a partir de la desviación estándar, indicando el número de veces que se la considera como medida de dispersión.

    Los puntos identificados en cada sector se seleccionan en 'grp': FIEC, FCNM, RECT.

    El tipo de medición tomada, 'tip', se identifica por: punto, 1m, gtw, dispositivo.

    Un parámetro auxiliar es 'LOS', que indica los puntos seleccionados con Línea de vista (1) y sin linea de vista (0). Para incluir todos de debe ingrear [1,0]. Este parametro se puede modificar en el archivo de entrada: arch_medidaubica.

    Los datos de cada eje se seleccionan mediante la función pares_usar(tabla, baliza, analiza, unabaliza, medida, modo) que entrega como resultado los arreglos de pares ordenados y las etiquetas con los nombres, par_etiqueta).

    La linealización se realiza con el método de los mínimos cuadrados, con lo que se establece el |error| promedio y desviación estándar.

    |error| = |yi - f(xi)| |error_{medio}| = \frac{1}{n}\sum|yi - f(xi)|

    Procedimiento aplicado

    Para el análisis primero se consideran todos los puntos disponibles para obtener la primera ecuación, mostrada en el ejemplo con la línea azul.

    Con ésto es posible determinar un error de estimación, para luego proceder a discriminar los puntos atípicos.

    Se realiza una nueva estimación de linealización habiendo discriminado los puntos atípicos y se observa el resultado.

    Resultados para baliza: gtwFIEC

    El resultado del algoritmo se presenta como gráfica, en pantalla y un archivo con los datos de las fórmulas.

    los resultados se pueden observar en lo mostrado.:

    baliza:  gtwFIEC
    Puntos usados: todos
    $ -10(4.908).log_{10}(d)+(1.406)$
    |error| promedio:  4.84  , std: 5.56
    Puntos usados: NoAtipico
    $ -10(5.12).log_{10}(d)+(6.714)$
    |error| promedio:  2.98  , std: 3.31
    >>> 
    

    Se observa que los valores fuera de la banda de valores con una desviación estándar (σ) se muestran distruidos en tres grupos: dos grupos a la izquierda y derecha de la gráfica por debajo de la banda y un grupo en el centro por sobre la banda.

    Se considera explorar la división del intervalo en dos, puesto que existen dos entornos: uno principalmente conformado con vegetación y otro con edificaciones.

    los resultados que se van al archivo, incluyen todos los decimales:

     exportar resultados :
    {'todos': {'intervalox': [52.543, 397.148], 
          'intervaloy': [-129.132183908046, -86.98969072164948], 
          'alpha': 4.907571379870146, 'beta': 1.4062384027235748, 
          'error_medio': 4.840153103044936, 'error_std': 5.562835152792785, 
          'eq_latex': '$ -10(4.908).log_{10}(d)+(1.406)$'
          }, 
      'NoAtipico': {'intervalox': [78.492, 397.148], 
          'intervaloy': [-129.132183908046, -86.98969072164948], 
          'alpha': 5.119532447831607, 'beta': 6.713572849706863, 
          'error_medio': 2.9780010745912833, 'error_std': 3.312804227070313, 
          'eq_latex': '$ -10(5.12).log_{10}(d)+(6.714)$'
          }
    } 
    

    Para revisar la situación se presentan los resultados con otra baliza.

    Baliza: gtwFCNM

    Resultados del algoritmo.

    baliza:  gtwFCNM
    Puntos usados: todos
    $ -10(5.403).log_{10}(d)+(8.423)$
    |error| promedio:  4.59  , std: 5.48
    Puntos usados: NoAtipico
    $ -10(5.574).log_{10}(d)+(10.758)$
    |error| promedio:  2.31  , std: 2.79
    

    Baliza: gtwRECT

    baliza:  gtwRECT
    Puntos usados: todos
    $ -10(4.89).log_{10}(d)+(8.541)$
    |error| promedio:  2.9  , std: 3.73
    Puntos usados: NoAtipico
    $ -10(4.587).log_{10}(d)+(0.326)$
    |error| promedio:  1.41  , std: 1.72
    

    Algoritmo en Python

    El algoritmo realiza el proceso de datos para cada baliza usando los datos del archivo "resumen_rssiUbica01.txt", que el el resultado del proceso realizado en Integrar las tablas de Rssi y coordenadas de los puntos

    Los resultados del algoritmo se almacenan en el archivo "arch_ecuaciones".

    Los parámetros para el análisis se incorporan en el diccionario "analiza". Los parámetros se describen al inicio de la página.

    Como el proceso de linealización se reutiliza, se lo incorpora como parte de la librería girni_lora_libreria, sin embargo la función se describe en detalle en Rssi(distancia) Linealización - función Python .

    Procedimiento

    Los datos se leen desde el archivo y se incorporan a una estructura de datos en Pandas.

    Para cada baliza se determina si se ha indicado 'analizar', con lo que se seleccionan los pares ordenados y etiquetas a usar mediante la función girni.pares_usar().

    Con los datos seleccionados, se aplica mínimos cuadrados  y se obtienen los errores mediante la función girni.linealiza_lstsq(). Mediante el criterio de desviación estándar se discriminan los datos atípicos y se vuelve a evaluar los datos sin atipicos, entregando el resultado mediante archivos y gráficas.

    # LoRa-Multipunto, Rssi vs distancia
    # linealización Rssi vs log10(distancia)
    # por mínimos cuadrados, Graficas 2D y 3D
    # Girni 2020-10-07 propuesta: edelros@espol.edu.ec
    
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    
    import girni_lora_libreria as girni
    
    # INGRESO
    # archivos de entrada
    modo = 'rx'
    medida = 'rssi'
    arch_medidaubica = 'resumen_rssiUbica01.txt'
    
    # archivos de salida
    arch_ecuaciones  = 'resumen_ecuacionSimple05.json'
    
    analiza = {'gtwRECT':{'analizar'  : 1,
                          'atipico_std' : 1,
                          'grp' : ['FIEC','RECT'],
                          'tip' : ['punto'],
                          'LOS' : [1,0]},
               'gtwFIEC':{'analizar'  : 1,
                          'atipico_std' : 1,
                          'grp' : ['FIEC','FCNM'],
                          'tip' : ['punto'],
                          'LOS' : [1,0]},
               'gtwFCNM':{'analizar'   : 1,
                          'atipico_std' : 1,
                          'grp' : ['FIEC','FCNM'],
                          'tip' : ['punto'],
                          'LOS' : [1,0]}
               }
    
    baliza = {'d1':'gtwRECT',
              'd2':'gtwFIEC',
              'd3':'gtwFCNM'}
    
    # Parámetros de grafica
    tipograf   = '2D'  # '2D','3D'
    escala     = 'log' # 'normal','log'
    escalabase = 10    # 10, np.exp()
    
    # PROCEDIMIENTO
    # Resultados de análisis
    ecuacion  = {}
    eq_graf = {}
    
    # leer datos
    tabla = pd.read_csv(arch_medidaubica, index_col='etiqueta')
    tabla = pd.DataFrame(tabla)
    
    # Analizar datos hacia una baliza
    for unabaliza in analiza:
    
        # Parámetros 
        analizar = analiza[unabaliza]['analizar']
        atipico_std = analiza[unabaliza]['atipico_std']
    
        if analizar:
            ecuacion[unabaliza] ={}
            eq_graf[unabaliza] = {}
            # pares a usar
            [pares,par_etiqueta] = girni.pares_usar(tabla,baliza,
                                                    analiza,unabaliza,
                                                    medida,modo)
            # analiza puntos para mínimos cuadrados
            xi = pares[:,0]
            yi = pares[:,1]
            
            ecuacion0 = girni.linealiza_lstsq(xi,yi)
    
            fdist0 = ecuacion0['eq_lambda']
            yi0  = fdist0(xi)
    
            # Selecciona atipicos
            dyi0std = ecuacion0['error_std']
            dyi0 = yi - yi0
            atipicos = np.abs(dyi0) >= dyi0std*atipico_std
            xi0_e = xi[atipicos]
            yi0_e = yi[atipicos]
            etiq0_e = par_etiqueta[atipicos]
    
            # datos sin atipicos ----------
            atipicoNo = np.abs(dyi0) <= dyi0std*atipico_std
            xi1 = xi[atipicoNo]
            yi1 = yi[atipicoNo]
            etiq1 = par_etiqueta[atipicoNo]
    
            ecuacion1 = girni.linealiza_lstsq(xi1,yi1)
    
            fdist1 = ecuacion1['eq_lambda']
            yi1  = fdist1(xi)
    
            # para exportar
            ecuacion[unabaliza] = {'todos': ecuacion0,
                                   'NoAtipico': ecuacion1
                                   }
            
            eq_graf[unabaliza]  = {'puntos': [xi,yi],
                                   'todos' : yi0,
                                   'atipicos':[xi0_e,yi0_e],
                                   'atip_etiq': etiq0_e,
                                   'NoAtipico':yi1
                                   }
    
    # SALIDA
    for unabaliza in ecuacion:
        print('baliza: ',unabaliza)
        for unaecuacion  in ecuacion[unabaliza]:
            error_medio = ecuacion[unabaliza][unaecuacion]['error_medio']
            error_std = ecuacion[unabaliza][unaecuacion]['error_std']
            print('Puntos usados:', unaecuacion)
            print(ecuacion[unabaliza][unaecuacion]['eq_latex'])
            print('|error| promedio: ',np.round(error_medio,2),
                  ' , std:',np.round(error_std,2))
    
        print('\n',ecuacion[unabaliza],'\n')
        print()
    
    # salida a archivo
    ecuacion = pd.DataFrame.from_dict(ecuacion)
    ecuacion.to_json(arch_ecuaciones)
    
    # GRAFICAR
    # Referencias para gráfica
    grupo   = ['FIEC' ,'FCNM'  ,'RECT','CIRC']
    colores = ['green','orange','grey','magenta']
    tipo    = ['punto','1m' ,'gtw','dispositivo']
    marcas  = [    'o','D'  ,'D'  ,'*' ]
    
    mostrargrpeti = ['FIEC','FCNM','RECT']
    mostrartipeti = ['1m','gtw']
    
    for unabaliza in ecuacion:
        figura,grafica = plt.subplots()
        if escala == 'log':
            grafica.set_xscale(escala,base=escalabase)
    
        # todos los puntos
        [xi, yi] = eq_graf[unabaliza]['puntos']
        grafica.scatter(xi,yi,marker='.')
        fdtxt = ecuacion[unabaliza]['todos']['eq_latex']
        
        # linea con todos los puntos
        yi0 = eq_graf[unabaliza]['todos']
        grafica.plot(xi,yi0,color='blue', label = fdtxt)
        
        [xi0_e,yi0_e] = eq_graf[unabaliza]['atipicos']
        etiq0_e = eq_graf[unabaliza]['atip_etiq']
        
        # cotas de error
        atipico_std = analiza[unabaliza]['atipico_std']
        dyi0std = ecuacion[unabaliza]['todos']['error_std']
        
        grafica.plot(xi,yi0 + dyi0std*atipico_std,
                        color='blue',linestyle='dotted')
        grafica.plot(xi,yi0 - dyi0std*atipico_std,
                        color='blue',linestyle='dotted')
        # atipicos
        grafica.scatter(xi0_e,yi0_e, color='red')
        # atipicos etiquetas
        m = len(xi0_e)
        for i in range(0,m,1):
            grafica.annotate(etiq0_e[i],
                            (xi0_e[i],yi0_e[i]),)
        
        # linea Sin Atipicos
        yi1 = eq_graf[unabaliza]['NoAtipico']
        fdtxt1 = ecuacion[unabaliza]['NoAtipico']['eq_latex']
        grafica.plot(xi,yi1, color='orange', label = fdtxt1)
    
        # etiquetas y títulos
        grafica.legend()
        grafica.set_ylabel(medida+'_'+modo)
        grafica.set_xlabel('distancia')
    
        untitulo = unabaliza+': '+medida+'_'+modo + ' vs distancia'
        grafica.set_title(untitulo)
        grafica.grid(True,linestyle='dotted',
                     axis='x', which='both')
        
        plt.show()
    
  • 3. Procesa datos. Modelo de pérdidas en propagación LoRa

    Para el modelo de perdidas de propagación, en cada punto se registra en el archivo las mediciones de Rssi y SNR. Cada archivo de datos procesan, tabulando y ordenanto los valores representativos del comportamiento del RSSI y SNR para revisar sus descriptores de estadística.

    Para realizar el procesamiento de los datos, se crearon algunas funciones y procedimientos para simplificar la escritura de instrucciones, las que se resumen en el archivo girni_lora_libreria.

    Las coordenadas geográficas de dada punto se registraron con un GPS diferencial usando el formato UTM en un archivo tipo texto.

    El procesamiento de los datos ser realiza en varios pasos donde se revisan los resultados parciales.

    El primero de ellos consiste en tabular los datos de Rssi y SNR de cada punto en un solo archivo, luego se añaden las coordenadas y distancias cada punto medido, para finalmente integrar ambos resultados en un solo archivo con RSSI, distancias, coordenadas de cada punto.

    Cada sección permite disponer de archivos intermedios que pueden ser usados para observar y procesar resultados que permitan realizar observaciones y mejoras a los modelos planteados. Entre los pasos intemedios está por ejemplo: observar en gráficas las ubicaciones de los puntos en el plano XY usando sus coordenadas, o en otro caso observar los valores de Rssi distribuidos en el espacio formado por el plano del ejemplo anterior y en el eje Z los valores promedios RSSI.

    Las siguientes secciones describen lo realizado para:

    Procesa datos. Resumen de archivos.txt por punto, Rssi o SNR

    Funciones girni_lora_libreria

    Procesa datos. Coordenadas GPS y distancias al archivo.txt

    Procesa datos. Ubica los puntos en gráfica 2D o 3D

    Procesa datos. Integra tablas Rssi, distancias y coordenadas

    Procesa Datos. Un punto - revisa descriptores estadísticos

  • 3.5 Procesa datos. Integra tablas Rssi, distancias y coordenadas

    Se integran los archivos de Rssi y coordenadas en un solo archivo. Este proceso permite procesar Rssi vs distancias para generar el modelo de la ecuación que los describe.

    Verificar los parámetros con los que se integran las tablas: nombres de archivos de entrada y salida, medidas.

    Ejemplo de resultados a obtener:

    registros:  57
    indices:  Index(['grupo', 'tipo', 'LOS_d1', 'LOS_d2', 'LOS_d3',
           'rssi_rx_d1','rssi_rx_d2', 'rssi_rx_d3', 'rssi_tx_d1', 'rssi_tx_d2', 'rssi_tx_d3',
           'c_norte', 'c_este', 'altitud', 'dist_d1', 'dist_d2', 'dist_d3','longitud', 'latitud'],
            dtype='object')
    
     Ejemplo de tabla: 
             grupo   tipo  LOS_d1  LOS_d2  LOS_d3  rssi_rx_d1  ...  altitud  dist_d1  dist_d2  dist_d3   longitud latitud
    etiqueta                                                   ...                                                         
    FIEC101   FIEC  punto       0       1       1 -123.096386  ...   62.213  423.450   78.492  351.924 -79.967500 -2.145463
    FIEC102   FIEC  punto       1       1       1 -120.094891  ...   82.308  451.423   67.435  357.220 -79.967762 -2.145404
    FIEC103   FIEC  punto       1       1       1 -116.923469  ...   77.421  449.414   60.141  364.706 -79.967697 -2.145337
    FIEC104   FIEC  punto       0       1       0 -116.046296  ...   87.352  399.437   97.611  343.328 -79.967307 -2.145563
    FIEC105   FIEC  punto       1       1       0 -119.713235  ...   82.708  382.250  117.785  326.571 -79.967231 -2.145730
    
    [5 rows x 19 columns]
    >>> 
    

    Archivo resumen de medida: resumen_rssimean01.txt

    archivo de resumen de coordenadas y distancias : resumen_ubica01.txt

    archivo de resultado integrado: resumen_RssiUbica01


    Algoritmo en Python

    La tabla general se constuye concatenando los componentes por colunnas.

    La primera parte se conforma con los datos del punto de las columnas de grupos, tipos y LOS que provienen del archivo de ubicaciones de los puntos.

    La segunda parte corresponde a la tabla de medidas de Rssi, y finalmente se complementa la tercera parte con la información de las coordenadas.

    La tabla resultante se almacena en un archivo de resumen.

    # Rssi y SNR LoRa punto a punto
    # LoRa-Multipunto, integra Ubicacion y Rssi
    # Girni 2020-10-07 propuesta: edelros@espol.edu.ec
    import numpy as np
    import pandas as pd
    
    # INGRESO
    # revisar parametros al inicio
    medida     = 'rssi'
    descriptor = 'mean'
    
    # archivos de entrada
    arch_rsmpuntos = 'resumen_rssimean02.txt'
    arch_rsmgps    = 'resumen_ubica01.txt'
    
    # archivos de salida
    arch_rssiubica = 'resumen_'+medida+'Ubica02.txt'
    
    # referencias
    baliza = {'d1':'gtwRECT',
              'd2':'gtwFIEC',
              'd3':'gtwFCNM'}
    
    # PROCEDIMIENTO
    # lectura de ubicacion
    ubica = pd.read_csv(arch_rsmgps, index_col='etiqueta')
    ubica = pd.DataFrame(ubica)
    n=len(ubica)
    
    # lectura de medida
    rsm_medida = pd.read_csv(arch_rsmpuntos)
    rsm_medida = pd.DataFrame(rsm_medida)
    rsm_medida.rename(columns={'Unnamed: 0':'etiqueta'},
                  inplace=True)
    rsm_medida['etiqueta'] = rsm_medida.etiqueta.astype(str)
    rsm_medida = rsm_medida.set_index('etiqueta')
    m=len(rsm_medida)
    
    # tabla concatenada por columnas Y UNION
    # para mantener referencias de coordenadas de vertices
    vertices  = list(baliza.keys())
    ultimovertice = vertices[-1]
    etiquetas = list(ubica.keys())
    donde = etiquetas.index('LOS_'+ultimovertice)
    
    # inicia con grupos, tipos y LOS
    tabla = ubica[etiquetas[:donde+1]]
    
    # continua con los datos de medida
    tabla = pd.concat([tabla,rsm_medida],
                      axis=1,join='outer')
    
    # completa con lo que resta de ubica
    tabla = pd.concat([tabla,
                       ubica[etiquetas[donde+1:]]],
                      axis=1,join='outer')
    tabla = tabla.rename_axis('etiqueta')
    k = len(tabla)
    
    # SALIDA
    print('registros: ',k)
    print('indices: ',tabla.keys())
    print('\n Ejemplo de tabla: ')
    print(tabla.head())
    tabla.to_csv(arch_rssiubica)
    
  • 3.4 Procesa datos. Ubica los puntos en gráfica 2D o 3D

    Para revisar los datos de coordenadas y ubicaciones obtenidos en el proceso anterior, se grafica la ubicación de los puntos usando las coordenadas este, norte y altura registradas con el GPS diferencial.

    Las gráficas se pueden realizar en 2D y en 3D, se adjunta los resultados:

    Como referencia se usan las etiquetas de los gateways que permiten observar la posición relativa de cada punto.

    Los grupos de puntos a mostrar se pueden seleccionar en el bloque de ingreso dentro de los parámetros de la gráfica, permitiendo obervar con mayor detalle la ubicación de cada punto

    Una observación en 3D de los puntos permite revisar que lso valores de altura presentan menor precisión que la de posición, principalmente en puntos ubicados entre vegetación.

    Se observa que las medidas de alturas tienen valores muy variables, por lo que el uso de los valores de altura se descarta inicialmente de los modelos. Se podría realizar posteriormente un análisis más detallado de lo presentado con los valores de las alturas.

    archivo de resumen de coordenadas y distancias : resumen_ubica01.txt


    Algoritmo en Python

    En el caso de gráficas 3D, para crear el archivo.gif animado es necesario instalar imagemagic. El nombre del archivo animado es 'rotando3D.gif' que se guarda en el mismo directorio del algoritmo.py.

    # Datos desde GPS Diferencial archivo.txt
    # Graficas de coordenadas 2D o 3D
    # Girni 2020-10-07 edelros@espol.edu.ec
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib.animation as animation
    
    # INGRESO
    # Archivo procesado con distancias
    arch_gpsrsm = 'resumen_ubica01.txt'
    
    # Parametros de gráfica
    mostrargrp = ['FIEC','FCNM','RECT']
    # ['FIEC','FCNM','RECT','CIRC']
    
    mostrartip = ['punto','1m','dispositivo','gtw']
    # ['punto','1m','dispositivo','gtw']
    
    tipograf = '3D' # '2D','3D'
    arch_rotacion = 'rotando3D.gif'
    
    # Referencias
    baliza  = {'d1':'gtwRECT',
               'd2':'gtwFIEC',
               'd3':'gtwFCNM'}
    grupo   = ['FIEC' ,'FCNM'  ,'RECT','CIRC']
    colores = ['green','orange','grey','magenta']
    tipo    = ['punto','1m' ,'gtw','dispositivo']
    marcas  = [    'o','D'  ,'D'  ,'*' ]
    
    # PROCEDIMIENTO
    # leer coordenadas
    ubica = pd.read_csv(arch_gpsrsm, index_col='etiqueta')
    ubica = pd.DataFrame(ubica)
    n = len(ubica)
    
    # puntos de vertices
    vertices = ubica.loc[baliza.values()]
    
    # segmento de grupo/tipo
    ubica['usar']  = False
    ubica['color'] = 'yellow'
    ubica['marca'] = 'o'
    for fila in ubica.index:
        unaeti  = str(fila)
        ungrupo = ubica['grupo'][fila]
        untipo  = ubica['tipo'][fila]
        cond1  = ungrupo in mostrargrp
        cond2  = untipo in mostrartip
        if cond1:
            cual = grupo.index(ungrupo)
            ubica.loc[fila,'color'] = colores[cual]
        if cond2:
            cual = tipo.index(untipo)
            ubica.loc[fila,'marca'] = marcas[cual]
        if (cond1 and cond2):
            ubica.loc[fila,'usar']  = True
    
    # SALIDA
    print(ubica.head())
    
    # Grafica 2D --------------------------
    if tipograf == '2D':
        figura, grafica = plt.subplots()
        for fila in ubica.index:
            unalon = ubica['c_este'][fila]
            unalat = ubica['c_norte'][fila]
            unaalt = ubica['altitud'][fila]
            usar   = ubica['usar'][fila]
            uncolor = ubica['color'][fila] 
            unamarca = ubica['marca'][fila]
            untipo = ubica['tipo'][fila]
            if usar:
                grafica.scatter(unalon,unalat,
                                color = uncolor,
                                marker = unamarca,
                                label = fila)
            if usar and untipo=='gtw':
                grafica.annotate(fila,
                                 (unalon,unalat))
        plt.xlabel('UTM_este')
        plt.ylabel('UTM_norte')
        plt.title('Ubicacion UTM')
        plt.show()
    
    # Grafica 3D --------------------------
    if tipograf == '3D':
        figura = plt.figure()
        grafica = Axes3D(figura)
        for fila in ubica.index:
            unalon  = ubica['c_este'][fila]
            unalat  = ubica['c_norte'][fila]
            unaalt  = ubica['altitud'][fila]
            usar    = ubica['usar'][fila]
            uncolor = ubica['color'][fila] 
            unamarca = ubica['marca'][fila]
            untipo  = ubica['tipo'][fila]
            if usar:
                grafica.scatter(unalon,unalat,
                                unaalt,
                                marker = unamarca,
                                color = uncolor,
                                label = fila)
            if (usar and (untipo=='gtw')):
                grafica.text(unalon,unalat,
                             unaalt,fila)
        grafica.set_xlabel('UTM_este')
        grafica.set_ylabel('UTM_norte')
        grafica.set_zlabel('altitud')
        grafica.set_title('Ubicacion UTM')
    
        def rotate(angle):
            grafica.view_init(azim=angle)
    
        print("realizando animation")
        rot_animation = animation.FuncAnimation(figura,
                        rotate,
                        frames = np.arange(45,360+45,10),
                        interval = 200)
        rot_animation.save(arch_rotacion, dpi=80,
                        writer = animation.PillowWriter(fps=5))
        plt.show()
    

    Referencia: Graficas 3D puntos dispersos-scatter en métodos numéricos

  • 3.3 Procesa datos. Coordenadas GPS y distancias al archivo.txt

    Procesa las coordenadas de cada punto medido, registradas en un archivo de texto usando un GPS Diferencial.

    Las coordenadas se encuentran en formato UTM ubicadas en la zona "17 M".

    Ejemplo de archivo de coordenadas del GPS Diferencial

    item,c_norte,c_este,altitud,etiqueta
    6,9762822.08,614817.64,62.213,FIEC101
    8,9762828.622,614788.477,82.308,FIEC102
    7,9762836.027,614795.73,77.421,FIEC103
    

    GPS Diferencial 01A los datos de cada punto se añaden los cálculos de distancia hacia cada "baliza". En el proceso, también se usa el nombre del archivo para separar el grupo y tipo de cada medición.

    Para observar los datos en el mapa mediante google Earth, se convierten las coordenadas UTM a latitud y longitud. Los valores de grupo permiten segmentar los puntos en el mapa.

    Referencias: https://en.wikipedia.org/wiki/Differential_GPS,
    https://es.wikipedia.org/wiki/Sistema_de_coordenadas_universal_transversal_de_Mercator


    Parámetros del algoritmo

    En el bloque de ingreso del algoritmo se configuran los parámetros que son requeridos para obtener un archivo más completo de coordenadas y distancias.

    Las posiciones para los gateways se indican como baliza, asociando el identificador d1, d2 o d3 con los nombres registrados con el GPS.

    En la tabla creada se usa la etiqueta de cada punto como índice de fila. las columnas corresponden a cada dato del punto.

    Un ejemplo de resultado al ejecutar del algoritmo es:

    archivo resumen: resumen_ubica01.txt
    una muestra de archivo: 
             grupo   tipo  LOS_d1  LOS_d2  LOS_d3      c_norte  ...  altitud  dist_d1  dist_d2  dist_d3  Longitud Latitud
    etiqueta                                                    ...                                                         
    FIEC101   FIEC  punto       1       1       1  9762822.080  ...   62.213  423.450   78.492  351.924 79.967500 -2.145463
    FIEC102   FIEC  punto       1       1       1  9762828.622  ...   82.308  451.423   67.435  357.220 79.967762 -2.145404
    FIEC103   FIEC  punto       1       1       1  9762836.027  ...   77.421  449.414   60.141  364.706 79.967697 -2.145337
    FIEC104   FIEC  punto       1       1       1  9762810.917  ...   87.352  399.437   97.611  343.328 79.967307 -2.145563
    FIEC105   FIEC  punto       1       1       1  9762792.518  ...   82.708  382.250  117.785  326.571 79.967231 -2.145730
    
    

    archivo de gps utm: ubicapuntos01.txt

    archivo de resumen de coordenadas y distancias : resumen_ubica01.txt


    Algoritmo en Python

    Para la conversión del sistema de coordenadas se usa la librería utm. En caso de no disponer de la librería, puede ser instalada con la instrucción pip desde una ventana de comandos:

    pip install utm

    con lo que es posible hacer las conversiones de UTM a latitud y longitud y viceversa.

    La fórmula de distancia en UTM es la tradicional para distancia entre dos puntos en el plano.

    # Datos desde GPS Diferencial a un archivo.txt
    # Incorpora distancias a los vertices (baliza),
    # añade grupo por sector, tipo de punto y LOS.
    # Girni 2020-10-07 edelros@espol.edu.ec
    import numpy as np
    import pandas as pd
    import utm
    
    # INGRESO
    # Archivo de coordenadas
    arch_gps    = 'ubicaPuntos01.txt'
    
    # Archivo salida procesado con distancias
    arch_gpsrsm = 'resumen_ubica03.txt'
    zona ='17 M'
    
    # Referencias
    baliza = {'d1':'gtwRECT',
              'd2':'gtwFIEC',
              'd3':'gtwFCNM'}
    grupo = ['FIEC','FCNM','RECT','CIRC']
    tipo  = ['punto','1m','gtw','dispositivo']
    
    # digitos decimales en distancias
    digitos = 3
    
    # PROCEDIMIENTO
    # leer coordenadas
    ubica = pd.read_csv(arch_gps,index_col='etiqueta')
    ubica.drop('item',inplace=True, axis=1)
    ubica = pd.DataFrame(ubica)
    n = len(ubica)
    zonanum = int(zona[0:2])
    zonalet = zona[3]
    
    # vertices con balizas
    baliza_key = list(baliza.keys())
    baliza_val = list(baliza.values())
    vertices   = ubica.loc[baliza_val]
    
    # distancias a vertices
    for fila in vertices.index:
        x1 = vertices['c_este'][fila]
        y1 = vertices['c_norte'][fila]
        x2 = ubica['c_este']
        y2 = ubica['c_norte']
        dist  = np.sqrt((x2-x1)**2 + (y2-y1)**2)
        dist  = np.round(dist,digitos)
        donde = baliza_val.index(fila)
        cual  = baliza_key[donde]
        etiq_col = 'dist_'+cual
        ubica[etiq_col] = dist
    
    # redondear a dos digitos
    ubica['c_este']  = np.round(ubica['c_este'],digitos)
    ubica['c_norte'] = np.round(ubica['c_norte'],digitos)
    ubica['altitud'] = np.round(ubica['altitud'],digitos)
    
    # añadir coordenadas en latitud y longitud
    ubica['longitud'] = 0.0
    ubica['latitud']  = 0.0
    for fila in ubica.index:
        lon_utm = ubica['c_este'][fila]
        lat_utm = ubica['c_norte'][fila]
        coord_gra = utm.to_latlon(lon_utm,lat_utm,
                                  zonanum,zonalet)
        ubica['latitud'][fila]  = coord_gra[0]
        ubica['longitud'][fila] = coord_gra[1]
    
    # grupo y tipo en cada punto
    ubica.insert(0,'grupo','')
    ubica.insert(1,'tipo','')
    for cadauno in ubica.index:
        
        # etiqueta de grupo
        esgrupo = 'CIRC'
        for ungrupo in grupo:
            cond1 = cadauno.startswith(ungrupo)
            cond2 = cadauno.endswith(ungrupo)
            if cond1 or cond2:
                esgrupo = ungrupo
        
        # etiqueta de tipo
        estipo = tipo[0]
        for untipo in tipo:
            cond1 = cadauno.startswith(untipo)
            if cond1:
                estipo = untipo
        if esgrupo == 'CIRC':
            estipo = 'dispositivo'
        
        ubica.loc[cadauno,'grupo'] = esgrupo
        ubica.loc[cadauno,'tipo']  = estipo
    
    # generar columna de Linea de vista LOS
    columna = 2
    for cadauno in baliza.keys():
        ubica.insert(columna,'LOS_'+cadauno,1)
        columna = columna + 1
    
    # SALIDA
    print('archivo resumen:', arch_gpsrsm)
    print('una muestra de archivo: ')
    print(ubica.head())
    print('vertices: \n',vertices)
    ubica.to_csv(arch_gpsrsm)
    
  • 3.2 Funciones girni_lora_libreria

    El archivo contiene un grupo de  funciones usadas para procesar los datos de este prototipo, se separaron del bloque principal de instrucciones con  el objetivo de simplificar el desarrollo del las actividades principales.

    Recuerde disponer de este archivo en la misma carpeta que el algoritmo que invoca a la librería.

    archivo de libreria: girni_lora_libreria

    datos de muestras

    tabulaPunto(unarchivo,directorio):
    Lee el archivo de un punto, tabula cada lectura por dispositivo en diccionario,
    por modo: rx, rx y remitente de cada paquete: baliza.

    describePunto(punto):
    estadistica descriptiva de un punto, calcula: count, mean, std, min, 25%, 50%, 75%, min, max

    resumen_medida(punto,medida,descriptor):
    Realiza la tabla resumen de puntos por la medida y descriptor a partir de la tabla de los puntos estadisticos descritos

    Linealización para ecuación

    pares_usar(tabla, baliza, analiza,unabaliza, unsector ='', medida = 'rssi',modo = 'rx')

    Selecciona desde la tabla puntos a usar respecto a una baliza, el resultado se entrega en una lista que contiene los pares ordenados y sus etiquetas [pares, par_etiqueta]

    linealiza_lstsq(xi,yi,digitos = 3)
    emplea el método de minimos cuadrados para entregar la ecuacion mediante un diccionario que contiene los parámetros para aplicarla.

    La variable dígitos indica cuántos dígitos se usarán para la ecuación en formato latex.

    El procedimiento se describe en: Rssi(distancia) Linealización – función Python

    unaecuacion = {'alpha'   : alpha,
                   'beta'    : beta,
                   'eq_latex': fdtxt0,
                   'intervalox' : [np.min(xi),np.max(xi)],
                   'error_medio': dyi0mean,
                   'error_std'  : dyi0std,
                   'eqg_latex'  : grtxt0,
                   'intervaloy'  : [np.min(yi),np.max(yi)],
                   'errorx_medio': dxi0mean,
                   'errorx_std'  : dxi0std,
                   }
    

    Proceso de triangulación

    dist_rssi(valor,ecuacion_rssi).  Evalua la ecuacion de distancia(rssi) revisando los intervalos disponibles para el valor de rssi dado.

    cruce2circulos(x1,y1,r1,x2,y2,r2). Revisa intervalo de area de cruce entre dos círculos de centro y radio: x1, y1, r1 // x2, y2, r2.

    raices2circulos(x1,y1,r1,x2,y2,r2,tolera=1e-10). Busca las intersección entre 2 circulos de centro y radio: x1,y1,r1 || x2,y2,r2 . Revisa con cruce2circulos().

    intersectacirculos(radio,centro,tolera=1e-10). Busca las intersecciones entre parejas de varios círculos y las entrega como [raicesx,raicesy], usa las funciones cruce2circulos() y con raices2circulos() que se encuentran en el enlace: Solución General de intersección de círculos

    trilatera(radio,centro,tolera=1e-10). Busca el baricentro entre las intersecciones de varios círculos, punto central en el área de
    intersección entre varios circulos. Requiere el resultado de la función:
    raiztodas = intersectacirculos(centro,radio, tolera = 1e-10)


    Algoritmo Python

    # Girni LoRa librerias 2020-10-07
    # LoRa-Multipunto, lecturas de Rssi y SNR
    # Girni 2020-10-07 propuesta: edelros@espol.edu.ec
    
    import numpy as np
    import pandas as pd
    import scipy.optimize as sp
    
    def tabulaPunto(unarchivo,carpeta, prefijo = 'multipunto'):
        ''' Lee el archivo de un punto dentro del carpeta,
            elimina el prefijo el nombre del archivo,
            tabula cada lectura por dispositivo en diccionario,
            por modo: rx, rx
            y remitente de cada paquete: baliza
            Prepara para procesar estadistica descriptiva
        '''
        # Datos estructura
        punto = {'rx':{},
                 'tx':{},
                 'nombre': ' '}
        
        # Lectura de unarchivo
        unarchivoubica = carpeta+'/'+unarchivo
        archivoPunto = open(unarchivoubica,'r')
    
        # nombre del punto desde nombre archivo
        pnombre = unarchivo
        pnombre = pnombre.strip('.txt')
        n = len(prefijo)
        pnombre = pnombre[n:]
        punto['nombre'] = pnombre
        
        linea = archivoPunto.readline()
        while (linea!=''):
            linea_rx = linea.startswith('rx')
            linea_tx = linea.startswith('tx')
            if linea_rx or linea_tx:
                # formato de trama:
                # tx_rx, c1_ff, d1_d2_d3, numtrama,
                # rssitx, snrtx, rssi_rx,snr_rx
                texto = linea.strip('\n')
                texto = texto.split(',')
                rx_tx      = texto[0]
                dir_recibe = texto[1]
                dir_remite = texto[2]
                ID_paquete = int(texto[3])
                rssi_tx    = float(texto[4])
                snr_tx     = float(texto[5])
                rssi_rx    = float(texto[6])
                snr_rx     = float(texto[7])
                # llena datos
                if dir_remite in punto[rx_tx].keys():
                    punto[rx_tx][dir_remite]['rssi_rx'].append(rssi_rx)
                    punto[rx_tx][dir_remite]['snr_rx'].append(snr_rx)
                    punto[rx_tx][dir_remite]['secuencia_rx'].append(ID_paquete)
                    punto[rx_tx][dir_remite]['rssi_tx'].append(rssi_tx)
                    punto[rx_tx][dir_remite]['snr_tx'].append(snr_tx)
                else:
                    punto[rx_tx][dir_remite] = {'rssi_rx': [rssi_rx],
                                                'snr_rx' : [snr_rx],
                                                'secuencia_rx':[ID_paquete],
                                                'rssi_tx': [rssi_tx],
                                                'snr_tx' : [snr_tx]}
            # siguiente línea
            linea = archivoPunto.readline()
        archivoPunto.close()
        return(punto)
    
    def describePunto(punto):
        ''' estadistica descriptiva de un punto
        calcula y registra count, mean, std,
        min,25%,50%,75%,max
        '''
        estadistica = {}
        for modo in punto.keys():
            estadistica[modo] = {}
            # analiza rssi y snr para rx y tx
            for disp in punto[modo].keys():
                estadistica[modo][disp] = ''
                valores = pd.DataFrame(punto[modo][disp])
                descrito = valores.describe()
                descrito = descrito.drop('secuencia_rx',axis=1)
                estadistica[modo][disp] = descrito       
        return(estadistica)
    
    
    def resumen_medida(punto,medida,descriptor):
        '''
        Realiza la tabla resumen de puntos por
        medida ('rssi' o 'snr')
        y descriptor ('count, mean, std,
         min,25%,50%,75%,max)
        a partir de la tabla de los puntos estadisticos descritos
        '''
        rsm_disp = pd.DataFrame()
        for disp in punto:
            undisp_rx = punto[disp][medida+'_rx']
            etiquetarx = medida+'_rx_'+disp
            rsm_disp[etiquetarx] = undisp_rx.copy()
        for disp in punto:
            undisp_tx = punto[disp][medida+'_tx']
            etiquetatx = medida+'_tx_'+disp
            rsm_disp[etiquetatx] = undisp_tx.copy()
        unafila = rsm_disp.loc[descriptor]
        return(unafila)
    
    def pares_usar(tabla,baliza, analiza,
                    unabaliza, unsector ='', 
                    medida = 'rssi', modo = 'rx'):
        ''' Selecciona en tabla los puntos a usar
            respecto a una baliza  y sector
            resultado en [pares, par_etiqueta]
        '''
        
        # balizas referencia para analizar
        baliza_key = list(baliza.keys())
        baliza_val = list(baliza.values())
    
        donde = baliza_val.index(unabaliza)
        cualbaliza = baliza_key[donde]
    
        # banderas de uso y atipicos
        bal_sec = cualbaliza+unsector
        tabla['usar_'+bal_sec] = 0
        tabla['atip_'+bal_sec] = 0
            
        # Parametros
        if unsector == '':
            analizarque = analiza[unabaliza]
        else:
            analizarque = analiza[unabaliza][unsector]
        atipico_std = analizarque['atipico_std']
        bal_grp = analizarque['grp']
        bal_tip = analizarque['tip']
        bal_LOS = analizarque['LOS']
                
        # usar segmento de grupo/tipo, bandera True/False
        for cadapunto in tabla.index:
            cond1 = tabla['grupo'][cadapunto] in bal_grp
            cond2 = tabla['tipo'][cadapunto] in bal_tip
            cond3 = tabla['LOS_'+cualbaliza][cadapunto] in bal_LOS
            cond4 = True
            if unsector != '':
                cond4 = tabla['sector_'+cualbaliza][cadapunto] == int(unsector.strip('s'))
            usar = cond1 and cond2 and cond3 and cond4
            valor = 0
            if usar:
                valor = 1
            tabla.loc[cadapunto,'usar_'+bal_sec] = valor
    
        # datos hacia baliza
        pares = []
        par_etiqueta = []
        for cadapunto in tabla.index:
            columna = medida+'_'+modo+'_'+cualbaliza
            xk = tabla['dist_'+cualbaliza][cadapunto]
            yk = tabla[columna][cadapunto]
            
            # no vacio y para usar
            cond1 = not(np.isnan(yk))
            cond2 = tabla['usar_'+bal_sec][cadapunto]
            if cond1 and cond2:
                unpar = np.array([xk,yk])
                
                # llena pares y etiquetas
                if len(pares)>0:
                    pares = np.concatenate((pares,[unpar]),axis=0)
                    par_etiqueta = np.concatenate((par_etiqueta,[cadapunto]),
                                                  axis=0)
                else:
                    pares = np.array([unpar])
                    par_etiqueta = np.array([cadapunto])
    
        # ordena pares para gráfica
        if len(pares)>0:
            ordenar = np.argsort(pares[:, 0])
            pares = pares[ordenar]
            par_etiqueta = par_etiqueta[ordenar]
    
        return ([pares, par_etiqueta])
    
    def linealiza_lstsq(xi,yi,digitos = 3):
        ''' usa minimos cuadrados para entregar la ecuacion
            digitos: usados en expresion latex
        '''
        unaecuacion = {}
        # Eje x en log10()
        xilog = np.log10(xi)
        n = len(xi)
        
        # mínimos cuadrados (least square),
        # distancia vs medida
        A = np.vstack([xilog, np.ones(n)]).T
        [m0, b0] = np.linalg.lstsq(A, yi, rcond=None)[0]
        alpha = -m0/10
        beta  = b0
    
        # ecuaciones expresion rssi(d)
        fdist0 = lambda d: -10*alpha*(np.log10(d))+beta
        
        fdtxt0 = r'$ rssi = -10(' + str(np.round(alpha,digitos))
        fdtxt0 = fdtxt0 + ')log_{10}(d)' # +('
        texto = '+'
        if beta <0:
            texto = '-'
        fdtxt0 = fdtxt0 + texto + str(np.round(np.abs(beta),digitos))+' $'
    
        # Errores respecto a rssi(d) 
        yi0  = fdist0(xi)
        dyi0 = yi - yi0
        dyi0mean = np.mean(np.abs(dyi0))
        dyi0std  = np.std(dyi0, dtype=np.float64)
    
        # ecuaciones expresion d(rssi)
        grssi0 = lambda rssi: 10**((beta-rssi)/(10*alpha))
        grtxt0 = r"$ d = 10^{(" + str(np.round(beta,digitos)) + ' - '
        grtxt0 = grtxt0 + 'rssi)/' + '(10('+str(np.round(alpha,digitos))+'))} $'
    
        # Errores respecto a rssi(d) 
        xi0  = grssi0(yi)
        dxi0 = xi - xi0
        dxi0mean = np.mean(np.abs(dxi0))
        dxi0std  = np.std(dxi0, dtype=np.float64)
        
        unaecuacion = {'alpha'   : alpha,
                       'beta'    : beta,
                       'eq_latex': fdtxt0,
                       'intervalox' : [np.min(xi),np.max(xi)],
                       'error_medio': dyi0mean,
                       'error_std'  : dyi0std,
                       'eqg_latex'  : grtxt0,
                       'intervaloy' : [np.min(yi),np.max(yi)],
                       'errorx_medio': dxi0mean,
                       'errorx_std'  : dxi0std,
                       }
        return(unaecuacion)
    
    def dist_rssi(valor,ecuacion_rssi, desplazar = 0):
        ''' evalua ecuacion de distancia(rssi)
            revisando los intervalos disponibles
        '''
        # resultados
        distancia = np.nan
        e_mean = np.nan
        e_1std = np.nan
        e_2std = np.nan
        
        # Revisa intervalos en ecuacion
        interv_fuera = 0
        tamano = len(ecuacion_rssi)
        eq_cual = list(ecuacion_rssi.keys())
    
        # revisa si hay ecuaciones
        if tamano >= 1:
            # en intervalo de todos los puntos?
            i_eq = 'r0'
            a = ecuacion_rssi[i_eq]['intervaloy'][0]
            b = ecuacion_rssi[i_eq]['intervaloy'][1]
            desplaza = 0; valor0 = valor
            if 'desplaza' in list(ecuacion_rssi[i_eq].keys()):
                desplaza = ecuacion_rssi[i_eq]['desplaza']
            if desplazar == 1:
                valor0 = valor - desplaza
            if valor0<a:
                interv_fuera = -1 # izquierda
            if valor0>=b:
                interv_fuera = 1  # derecha
        if interv_fuera!=0:
            if interv_fuera == 1: # derecha
                i_eq = eq_cual[1]
            if interv_fuera == -1: # izquierda
                i_eq = eq_cual[-1]
            # intervalo rssi [a,b)
            a = ecuacion_rssi[i_eq]['intervaloy'][0]
            b = ecuacion_rssi[i_eq]['intervaloy'][1]
            alpha = ecuacion_rssi[i_eq]['alpha']
            beta = ecuacion_rssi[i_eq]['beta']
            # correccion de formula por desplazamiento
            desplaza = 0; valor1 = valor
            if 'desplaza' in list(ecuacion_rssi[i_eq].keys()):
                desplaza = ecuacion_rssi[i_eq]['desplaza']
            if desplazar == 1:
                valor1 = valor - desplaza
            distancia = 10**((valor1-beta)/(-10*alpha))
            
            # errores estimados de distancia
            error_medio = ecuacion_rssi[i_eq]['error_medio']
            error_std  = ecuacion_rssi[i_eq]['error_std']
            dist_mean = 10**((valor-error_medio-beta)/(-10*alpha))
            dist_1std = 10**((valor-error_std-beta)/(-10*alpha))
            dist_2std = 10**((valor-2*error_std-beta)/(-10*alpha))
            e_mean = np.abs(distancia-dist_mean)
            e_1std = np.abs(distancia-dist_1std)
            e_2std = np.abs(distancia-dist_2std)
    
        # evalua valor si hay ecuacion
        if tamano>1 and interv_fuera == 0:
            donde = eq_cual.index('r0')
            eq_cual.pop(donde)
            # Revisa que exista ecuacion en cada intervalo
            for i_eq in eq_cual:
                if ecuacion_rssi[i_eq] is None:
                    donde = eq_cual.index(i_eq)
                    eq_cual.pop(donde)
                    
            # revisa intervalo y evalua
            for i_eq in eq_cual:
                # correccion de formula por desplazamiento
                desplaza = 0; valor1 = valor
                if 'desplaza' in list(ecuacion_rssi[i_eq].keys()):
                    desplaza = ecuacion_rssi[i_eq]['desplaza']
                if desplazar == 1:
                    valor1 = valor - desplaza
                # intervalo rssi [a,b)
                a = ecuacion_rssi[i_eq]['intervaloy'][0]
                b = ecuacion_rssi[i_eq]['intervaloy'][1]
                
                cond1 = (valor1>=a) and (valor1<b)
                
                if cond1:
                    alpha = ecuacion_rssi[i_eq]['alpha']
                    beta = ecuacion_rssi[i_eq]['beta']
                    distancia = 10**((valor1-beta)/(-10*alpha))
                    
                    # errores estimados de distancia
                    error_medio = ecuacion_rssi[i_eq]['error_medio']
                    error_std  = ecuacion_rssi[i_eq]['error_std']
                    dist_mean = 10**((valor-error_medio-beta)/(-10*alpha))
                    dist_1std = 10**((valor-error_std-beta)/(-10*alpha))
                    dist_2std = 10**((valor-2*error_std-beta)/(-10*alpha))
                    e_mean = np.abs(distancia-dist_mean)
                    e_1std = np.abs(distancia-dist_1std)
                    e_2std = np.abs(distancia-dist_2std)
        
        return([distancia,e_mean,e_1std,e_2std,interv_fuera])
    
    # Las siguientes funciones tienen como objetivo
    # realizar la trilateración entre varios circulos
    # def trilatera(radio,centro,tolera=1e-10)
    # requieren las funiones:
    # def cruce2circulos(x1,y1,r1,x2,y2,r2)
    # def raices2circulos(x1,y1,r1,x2,y2,r2,tolera=1e-10)
    # def intersectacirculos(radio,centro,tolera=1e-10)
    
    def cruce2circulos(x1,y1,r1,x2,y2,r2):
        ''' Revisa intervalo de area de cruce
            entre dos círculos de centro y radio
            x1,y1,r1 // x2,y2,r2
        '''
        uncruce = []
        dx = x2 - x1
        dy = y2 - y1
        d_centros = np.sqrt(dx**2 + dy**2)
        d_cruce   = r2 + r1
        
        # los circulos se cruzan o tocan
        if d_cruce >= d_centros:
    
            # intervalos de cruce
            xa = np.max([x1-r1,x2-r2])
            xb = np.min([x1+r1,x2+r2])
            ya = np.max([y1-r1,y2-r2])
            yb = np.min([y1+r1,y2+r2])
            
            # cada circulo arriba, abajo
            abajo1 = 0 ; arriba1 = 0
            abajo2 = 0 ; arriba2 = 0
            if ya<=y1:
                abajo1  = 1
            if yb>=y1:
                arriba1 = 1
            if ya<=y2:
                abajo2  = 1
            if yb>=y2:
                arriba2 = 1
            sector  = [ abajo1*abajo2, abajo1*arriba2,
                       arriba1*abajo2, arriba1*arriba2]
            uncruce = [xa,xb,ya,yb,sector]
        return(uncruce)
    
    def raices2circulos(x1,y1,r1,x2,y2,r2,tolera=1e-10):
        ''' busca las intersección entre 2 circulos
            de centro y radio: x1,y1,r1 || x2,y2,r2
            revisa con cruce2circulos()
        '''
        casicero = tolera*np.min([r1,r2])
        uncruce = cruce2circulos(x1,y1,r1,x2,y2,r2)
        raizx = []; raizy = []
        secruzan =  0
        
        # si hay cruce de circulos
        if len(uncruce)>0:
            sectores = [[-1,-1],[-1,1], 
                        [ 1,-1],[ 1,1]]
            [xa,xb,ya,yb,sector] = uncruce
            xc = (xa+xb)/2
            dx = xb-xa
            dy = yb-ya
            k = len(sector)
            if dx<casicero: # se tocan en un punto
                k = 1
            for j in range(0,k,1):
                if sector[j]==1:
                    s1 = sectores[j][0]
                    s2 = sectores[j][1]
                    def gx(x,x1,r1,casicero):
                        z = r1**2-(x-x1)**2
                        if np.abs(z)<casicero:
                            z = 0
                        return(z)
                    fx1 = lambda x: s1*np.sqrt(gx(x,x1,r1,casicero)) + y1
                    fx2 = lambda x: s2*np.sqrt(gx(x,x2,r2,casicero)) + y2
                    fx  = lambda x: fx1(x)-fx2(x)
                    
                    fa = fx(xa)
                    fb = fx(xb)
                    raiz1 = np.nan
                    raiz2 = np.nan
                    
                    # intervalo/2 izquierda
                    xc = xc + dx*tolera
                    fc = fx(xc)
                    cambio = np.sign(fa)*np.sign(fc)
                    if cambio<0:
                        raiz1 = sp.bisect(fx,xa,xc,xtol=tolera)
                        
                    # intervalo/2 derecha
                    xc = xc - 2*dx*tolera
                    fc = fx(xc)
                    cambio = np.sign(fc)*np.sign(fb)
                    if cambio<0:
                        raiz2 = sp.bisect(fx,xc,xb,xtol=tolera)
                        
                    # si hay contacto en un borde
                    if dx<casicero and dy>0:
                        raiz1 = xa
                    if dy<casicero and dx>0:
                        raiz1 = x1
                        
                    # Añade si existe raiz
                    if not(np.isnan(raiz1)):
                        raizx.append(raiz1)
                        raizy.append(fx1(raiz1))
                    if not(np.isnan(raiz2)):
                        raizx.append(raiz2)
                        raizy.append(fx1(raiz2))
                    secruzan = 1
        # No hay cruce de circulos
        if len(uncruce) == 0:
            dx = x2 - x1
            dy = y2 - y1
            m = dy/dx
            theta = np.arctan2(dy,dx)
            dx1 = r1* np.cos(theta)
            dx2 = r2* np.cos(theta)
            xi1 = x1 + dx1
            xi2 = x2 - dx2
            b = y1 - m*x1
            raizx = [(xi1+xi2)/2]
            raizy = [m*raizx[0]+b]
            
        raices = [raizx,raizy,secruzan]
        return(raices)
    
    def intersectacirculos(radio,centro,tolera=1e-10):
        ''' busca las intersecciones entre parejas de varios
            círculos y las entrega como [raicesx,raicesy]
            usa las funciones cruce2circulos()
            y con raices2circulos() que se encuentran en el enlace:
        http://blog.espol.edu.ec/analisisnumerico/s3eva_it2018_t1-interseccion-de-dos-circulos/#intersecta2circulos
        '''
        vertices = list(centro.keys())
        n = len (vertices)
        
        # agrupa raices en todasx y todasy
        todasx  = [] ; todasy  = []
        cruces = np.zeros(shape=(n,n),dtype=int)
        for i in range(0,n-1,1):
            for j in range(i+1,n,1):
                x1 = centro[vertices[i]][0]
                y1 = centro[vertices[i]][1]
                r1 = radio[vertices[i]]
                x2 = centro[vertices[j]][0]
                y2 = centro[vertices[j]][1]
                r2 = radio[vertices[j]]
    
                # busca raices entre 2 circulos
                raices = raices2circulos(x1,y1,r1,x2,y2,r2,tolera)
                raizx = raices[0]
                raizy = raices[1]
                cruces[i,j] = raices[2]
                cruces[j,i] = raices[2]
                m = len(raizx)
                if m>0:
                    for k in range(0,m,1):
                        todasx.append(raizx[k])
                        todasy.append(raizy[k])
                    
        raiztodas = [todasx,todasy,cruces]
        return(raiztodas)
    
    def trilatera(radio,centro,tolera=1e-10):
        ''' busca el baricentro entre las intersecciones
            de varios círculos punto central en el área de
            intersección entre varios circulos.
            requiere el resultado de la función:
            raiztodas = intersectacirculos(centro,radio, tolera = 1e-10)
        '''
        vertices = list(centro.keys())
        n = len (vertices)
        # revisa raiz dentro de cada circulo
        raiztodas = intersectacirculos(radio,centro,tolera)
        todasx = raiztodas[0]
        todasy = raiztodas[1]
        cruces = raiztodas[2]
        m = len(todasx)
        raicesx = []
        raicesy = []
        fuera = []
        for k in range(0,m,1):
            xk = todasx[k]
            yk = todasy[k]
            dentro = 0
            for i in range(0,n,1):
                x1 = centro[vertices[i]][0]
                y1 = centro[vertices[i]][1]
                r1 = radio[vertices[i]]
                dx = x1-xk
                dy = y1-yk
                d_centro = np.sqrt(dx**2+dy**2)
                if d_centro<=(r1*(1+tolera)):
                    dentro = dentro + 1
            if dentro == n:
                raicesx.append(xk)
                raicesy.append(yk)
                
        # busca baricentro
        baricentro = np.nan
        barerror = np.nan
        q = len(raicesx)
        if q>0:
            xbar = np.mean(raicesx)
            ybar = np.mean(raicesy)
            baricentro = [xbar,ybar]
            barerror = 0
            for i in range(0,q,1):
                d = np.sqrt((xbar-raicesx[i])**2+(ybar-raicesy[i])**2)
                if d>barerror:
                    barerror = d
            poligono = [raicesx,raicesy]
        else:
            poligono = [todasx,todasy]
        resultado = {'baricentro': baricentro,
                     'barerror'  : barerror,
                     'poligono'  : poligono,
                     'nocruzaen' : ''}
    
        
        # revisa espacio entre circulos
        sumacruces = np.sum(cruces,axis = 0)
        if 0 in list(sumacruces):
            raicesx = todasx.copy()
            raicesy = todasy.copy()
            for i in range(0,n,1):
                if sumacruces[i]==0:
                    x1 = centro[vertices[i]][0]
                    y1 = centro[vertices[i]][1]
                    lejanamax = -1
                    lejana_en = -1
                    for k in range(0,len(raicesx),1):
                        dx = x1 - raicesx[k]
                        dy = y1 - raicesy[k]
                        d_centro = np.sqrt(dx**2+dy**2)
                        if d_centro>lejanamax:
                            lejanamax = d_centro
                            lejana_en = k
                    if lejana_en>=0:
                        raicesx.pop(lejana_en)
                        raicesy.pop(lejana_en)
            # busca baricentro
            baricentro = np.nan
            barerror = np.nan
            q = len(raicesx)
            if q>0:
                xbar = np.mean(raicesx)
                ybar = np.mean(raicesy)
                baricentro = [xbar,ybar]
                barerror = 0
                for i in range(0,q,1):
                    d = np.sqrt((xbar-raicesx[i])**2+(ybar-raicesy[i])**2)
                    if d>barerror:
                        barerror = d
                poligono = [raicesx,raicesy]
            else:
                poligono = [todasx,todasy]
            resultado = {'baricentro': baricentro,
                         'barerror'  : barerror,
                         'poligono'  : poligono,
                         'nocruzaen' : vertices[i] }
    
    
        return(resultado)
    
  • 3.1 Procesa datos. Resumen de archivos.txt por punto, Rssi o SNR

    Los archivos de las capturas de datos por USB-Serial que se obtienen en las campañas de medición en cada punto, se requieren: procesar, ordenar y obtener los valores representativos del comportamiento del RSSI y SNR.

    Cada archivo contiene al menos 100 lecturas a cada baliza. Los archivos se encuentran en una carpeta desde donde se procesarán y se obtendrá un archivo "resumen".

    Para procesar los archivos de deben indicar los siguientes parámetros:

    • medida a observar: rssi, snr
    • descriptor estadistico a usar: mean, std, count, max, min, ...
    • trama: rx
    • una 'carpeta' donde se encuentran los archivos a procesar
    • nombre del archivo de resultados

    Un ejemplo de los resultados en el resumen por punto del promedio Rssi es:

    resumen trama: rx  medida: rssi
             rssi_rx_d1  rssi_rx_d2  rssi_rx_d3  rssi_tx_d1  rssi_tx_d2  rssi_tx_d3
    FCNM101 -124.396907 -129.132184  -81.620000  -98.000000 -126.632184  -78.040000
    FCNM102 -114.345455 -118.917526  -92.053571  -91.000000 -118.474227  -80.482143
    FCNM103 -115.459854 -124.634615  -95.453901  -91.000000 -123.173077  -76.617021
    FCNM104 -115.922414 -123.678899 -110.090909  -91.000000 -122.844037 -107.933884
    FCNM105 -108.239726 -124.893204 -102.007042  -91.000000 -123.582524  -99.281690
    FCNM106 -128.383929 -127.910000 -117.589041  -91.000000 -127.290000 -110.082192
    FCNM107 -124.429907 -123.009009 -105.078125  -91.000000 -123.225225 -102.343750
    

    Carpeta con archivos a procesar: lectura_puntos

    Archivos resultado:

    para Rssi: resumen_rssimean01.txt

    Para SNR: resumen_snrmean01.txt


    Algortimo en Python

    Para simplificar el algoritmo, se realizó una libreria de funciones, "girni_lora_libreria.py",  que se debe de ubicar en la misma carpeta del algoritmo que lo usa.

    Para almacenar los datos tabulados, se crea tambien el archivo de ésta tabla en formato .json, en caso de que se requiera realizar otro tipo de análisis a los datos.

    '''
    Procesa Archivo.txt de LoRa Rssi y SNR
    de muestras desde un puerto Serial
    Girni 2020-10-07 propuesta: edelros@espol.edu.ec
    medida: 'rssi', 'snr'
    unestadistico: 'count', 'mean', 'std',
          'min', '25%', '50%', '75%', 'max'
    modo: rx baliza a dispositivo
          tx difusion hacia balizas
    '''
    
    import os
    import numpy as np
    import pandas as pd
    import girni_lora_libreria as girni
    
    # INGRESO
    # revisar parametros al inicio
    medida     = 'rssi'
    descriptor = 'mean'
    trama      = 'rx'
    
    # carpeta de Archivos entrada
    carpeta = 'Lecturas_dispositivo'
    
    # Archivos para resultados
    arch_rsmUnestadistico = 'resumen_'+medida+descriptor+'02.txt'
    arch_detalle = 'tabla_puntosdatos2.json'
    
    # PROCEDIMIENTO
    # lista los archivos entrada en la carpeta
    archivosPuntos = os.listdir(carpeta)
    
    # Lee cada archivo de un punto y tabula
    # las lecturas por cada trama: rx, tx
    # cada remitente de paquete: baliza
    tabula = {}
    for unarchivo in archivosPuntos:
        unpunto  = girni.tabulaPunto(unarchivo,carpeta)
        p_nombre = unpunto['nombre']
        tabula[p_nombre] = unpunto
    tabula = pd.DataFrame(tabula)
    tabula = tabula.drop(['nombre'])
    tabula = tabula.T
    
    # Estadistica descriptiva de cada punto tabulado
    descrito = {}
    for cadapunto in tabula.index:
        punto = tabula.loc[cadapunto]
        descrito[cadapunto] = girni.describePunto(punto)
    
    # resumen de una medida: rssi, snr
    rsm_medida = pd.DataFrame()
    for cadapunto in descrito.keys():
        punto = descrito[cadapunto][trama]
        unafila = girni.resumen_medida(punto,medida,descriptor)
        rsm_medida[cadapunto] = unafila
    rsm_medida = rsm_medida.T
    rsm_medida = rsm_medida.sort_index(axis=1)
    
    # SALIDA
    # muestra en pantalla solo un ejemplo.head()
    print('resumen trama: ', trama,', medida: ', medida)
    print(rsm_medida.head(n=10))
    
    # Se escribe todo el archivo a csv
    rsm_medida.to_csv(arch_rsmUnestadistico)
    
    # Para grabar en formato .json
    tabula.to_json(arch_detalle)
    
  • 2.2 Captura datos. Dispositivo con USB-Serial a un archivo.txt usando Python

    La relación de valores de Rssi y distancia hacia cada punto de referencia (baliza), se estima a partir de las mediciones en varios puntos establecidos y ubicados en un mapa.

    En cada punto, los valores se estiman considerando usar al menos 100 lecturas hacia cada baliza en cada punto.

    En éste proceso usa un computador portátil que captura las tramas de cada paquete recibido por el dispositivo desde el puerto USB-Serial a un archivo.txt.

    Cada baliza periodicamente emite un paquete que permite al dispositivo registar los valores de Rssi y SNR. Cada captura se registra en una línea de texto como un paquete recibido y etiquetado como "rx".

    Para registrar el Rssi y SNR desde el otro extremo, la baliza, se emite un paquete de difusión desde el dispositivo, con el que la baliza al recibirlo toma los valores de Rssi y SNR y los envia al dispositivo en el próximo envío. El paquete de difusión del dispositivo se etiqueta como 'tx' en el archivo txt.

    En la pantalla del computador portátil se muesta el estado del registro de datos en cada evento de recepción.

    Configuración del algoritmo en Python

    La captura de datos desde USB se realiza usando instrucciones en Python, con la que se pre-procesa cada paquete y se almacena en un archivo.

    En el bloque de INGRESO es necesario indicar:

    • el puerto 'com' y los baudios del dispositivo conectado por USB
    • el nombre del archivo de texto que registra la captura de datos
    • las direcciones de las balizas, identificadores indicados en la sección anterior
    • la cantidad mínima de lecturas desde cada baliza

    El puerto 'com' se puede verificar antes de iniciar la captura, usando el IDE de arduino, al conectar lo muestra en el menu de herramientas/puerto.

    El proceso de captura se ejecuta en un lazo, en caso de ser necesario deterner se pueden usar las teclas "Ctrl-c".

    En cada punto de referencia se requiere asignar un nombre a cada archivo, iniciando el nombre con "multipunto" añadiendo letras para identificar el sector y un número como identificador del punto.

    En los puntos de prueba el nombre del archivo inicia con "mpcircuito". Observe que ambas palabras de inicio contienen la misma cantidad de caracteres, pues facilita el procesamiento de los archivos.

    Los archivos se almacenan en una carpeta de donde se puedan leer por lotes en el siguiente proceso de los datos.


    Algoritmo en Python

    ''' Rssi y SNR LoRa punto a punto
    paquete de datos desde puerto Serial
    generados desde dispositivo LoRa
    modelo Heltec Lora 32 v.2
    Girni 2020-10-07 propuesta: edelros@espol.edu.ec
    '''
    import numpy as np
    import matplotlib.pyplot as plt
    import serial, time
    
    # INGRESO
    # Puerto de captura de datos USB-Serial
    puerto = 'com3'
    baudios = 115200
    
    # Archivo para el registro de cada evento rx,tx
    nombrearchivo = 'multipuntoFIEC101.txt'
    
    # identificadores de balizas
    baliza = ['d1','d2','d3']
    
    # Mínimo de lecturas por baliza
    lecturasMin = 100
    
    # PROCEDIMIENTO
    # Crea registro de lecturas
    n = len (baliza)
    registro = {}
    for elemento in baliza:
        registro[elemento] = {'cuenta': 0 ,
                              'rssi': np.array([0.0,0.0]),
                              'snr':  np.array([0.0,0.0]),
                              'SumaRssi': np.array([0.0,0.0]),
                              'SumaSnr':  np.array([0.0,0.0]),
                              'minmaxRssi': np.array([0.0,-200.0,0.0,-200.0]),
                              'minmaxSnr': np.array([100.0,0.0,100.0,0.0])
                              }
    
    # inicializa archivo.txt a vacio
    archivo = open(nombrearchivo,'w')
    archivo.close()  # Cierra el archivo
    
    # Abre puerto Serial
    arduino = serial.Serial(puerto, baudios)
    arduino.setDTR(False)
    time.sleep(0.3)
    
    # limpia buffer de datos anteriores
    arduino.flushInput()  
    arduino.setDTR()  
    time.sleep(0.3)
    print('\nEstado del puerto: ',arduino.isOpen())
    print('Nombre del dispositivo conectado: ', arduino.name)
    print('Dump de la configuración:\n ',arduino)
    print('\n###############################################\n')
    
    # Lectura de datos
    np.set_printoptions(precision=2)
    conteo = 0
    difunde = 0
    while conteo<lecturasMin:
        #esperar hasta recibir un paquete
        while (arduino.inWaiting()==0):
            pass 
    
        # leer linea desde puerto serial
        linea = arduino.readline()
        # binario a texto, elimina /r/n
        texto = linea.decode()
        texto = linea.strip()
        
        # identificar la trama como rx, tx
        cond1 = texto.startswith('tx')
        cond2 = texto.startswith('rx')
        if cond1 or cond2:
            archivo = open(nombrearchivo,'a')
            archivo.write(texto+'\n')
            archivo.close()
            
            if (texto.startswith('tx')):
                difunde = difunde + 1
            if (texto.startswith('rx')):
                texto = texto.split(',')
                tipo  = texto[0]
                dir_remite = texto[2]
                paqrcbvID  = texto[3]
                rssi_tx    = float(texto[4])
                snr_tx     = float(texto[5])
                rssi_rx    = float(texto[6])
                snr_rx     = float(texto[7])
                if tipo == "rx":
                    # conteo de lecturas
                    cual = dir_remite
                    registro[cual]['cuenta']=registro[cual]['cuenta']+1
                    if registro[cual]['cuenta']<conteo:
                        conteo = registro[cual]['cuenta']
    
                    # acumulado
                    registro[cual]['SumaRssi'][0] = registro[cual]['SumaRssi'][0]+rssi_tx
                    registro[cual]['SumaSnr'][0]  = registro[cual]['SumaSnr'][0]+snr_tx
                    registro[cual]['SumaRssi'][1] = registro[cual]['SumaRssi'][1]+rssi_rx
                    registro[cual]['SumaSnr'][1]  = registro[cual]['SumaSnr'][1]+snr_rx
    
                    # promedios
                    cuantos = registro[cual]['cuenta']
                    registro[cual]['rssi'] = registro[cual]['SumaRssi']/cuantos
                    registro[cual]['snr']  = registro[cual]['SumaSnr']/cuantos
    
                    # minimos y maximos
                    registro[cual]['minmaxRssi'][0] = np.min([rssi_tx,registro[cual]['minmaxRssi'][0]])
                    registro[cual]['minmaxRssi'][1] = np.max([rssi_tx,registro[cual]['minmaxRssi'][1]])
                    registro[cual]['minmaxRssi'][2] = np.min([rssi_rx,registro[cual]['minmaxRssi'][2]])
                    registro[cual]['minmaxRssi'][3] = np.max([rssi_rx,registro[cual]['minmaxRssi'][3]])
    
                    registro[cual]['minmaxSnr'][0] = np.min([snr_tx,registro[cual]['minmaxSnr'][0]])
                    registro[cual]['minmaxSnr'][1] = np.max([snr_tx,registro[cual]['minmaxSnr'][1]])
                    registro[cual]['minmaxSnr'][2] = np.min([snr_rx,registro[cual]['minmaxSnr'][2]])
                    registro[cual]['minmaxSnr'][3] = np.max([snr_rx,registro[cual]['minmaxSnr'][3]])
    
                # Muestra en pantalla el estado de recepción
                print('\n difusion: ',difunde)
                print(texto)
                for elemento in baliza:
                    print(elemento,registro[elemento]['cuenta'],
                          '\tprom Rssi[tx,rx] \t   Snr[tx,rx]')
                    print("prom   :",registro[elemento]['rssi'],"\t  ",
                          registro[elemento]['snr'])
                    print("min,max:",registro[elemento]['minmaxRssi'],
                          registro[elemento]['minmaxSnr'])
    
    # Cerrar el puerto serial.
    serial.Serial.close