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()

4.4 Rssi vs Distancia. Linealiza POR intervalos – Algoritmo Python

El algoritmo para realizar dos o más intervalos considera usar un estimado de frontera a lo largo de las mediciones.

Se implementa añadiendo al bloque de ingreso dos parámetros en el bloque ‘analiza’:

1. frontera es un vector donde se indican las distancia de los puntos de corte de los subintervalos sin considerar los extremos, corresponden a la frontera estimada en el mapa. El algoritmo actualiza la frontera con la intersección  de la linealización de dos subintervalos consecutivos. Si frontera es un vector vacío [] se asume que se trabaja con todo el intervalo.

2. ‘atipInterv_std’ es un vector para la discriminación de los valores atipicos aplicada en cada subintervalo.

Como referencia para comparar, el resultado del análisis de todo el intervalo se denomina ‘r0’. Los subintervalos se identifican por ‘r1′,’r2’, etc, en orden al alejarse de la baliza.

En el resultado como valor complementario de revisión se añade el coeficiente de correlación de los puntos usados en el análisis.

Algoritmo en Python

# 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 json
import matplotlib.pyplot as plt
import girni_lora_libreria as girni

# INGRESO
# archivos de entrada
modo       = 'rx'
medida     = 'rssi'
descriptor = 'mean'
arch_medidaubica = 'rsmP06_'+medida+'Ubica01Intervalo.txt'

# archivos de salida
arch_ecuaciones  = 'rsmP07_ecuacion01Intervalos.json'

# Analizar por segmentos
analiza = {'gtwRECT':{'analizar'  : 1,
                      'atipico_std' : 1,
                      'frontera'    : [320],
                      'atipInterv_std': [1,1],
                      'p_amplia': 4,
                      'grp' : ['RECT','FIEC'],
                      'tip' : ['punto'],
                      'LOS' : [1]},
           'gtwFIEC':{'analizar'  : 1,
                      'atipico_std' : 1,
                      'frontera'    : [190], 
                      'atipInterv_std': [1,1],
                      'p_amplia': 4,
                      'grp' : ['FIEC','FCNM'],
                      'tip' : ['punto'],
                      'LOS' : [1,0]},
           'gtwFCNM':{'analizar'   : 1,
                      'atipico_std' : 1,
                      'frontera'   : [235.0],
                      'atipInterv_std': [2,2],
                      'p_amplia': 4,
                      '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()
casicero   = 1e-4
precision  = 3
intersectar = 1 # 0:Falso, 1: Verdadero

# Referencias de 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']

# PROCEDIMIENTO
# Resultados de análisis
ecuacion = {}
eq_graf  = {}

# leer datos
tabla = pd.read_csv(arch_medidaubica, index_col='etiqueta')
tabla = pd.DataFrame(tabla)

# analiza datos hacia una baliza
for unabaliza in analiza:

    # Parámetros
    analizar = analiza[unabaliza]['analizar']

    if analizar:
        ecuacion[unabaliza] = {}
        eq_graf[unabaliza]  = {}
        
        # pares a usar de baliza
        [pares,par_etiqueta] = girni.pares_usar(tabla,baliza,
                                            analiza,unabaliza,
                                            medida = medida ,
                                            modo = modo)
        # todos los puntos
        # analiza puntos para mínimos cuadrados
        xi = pares[:,0]
        yi = pares[:,1]
        n_xi = len(xi)
        # coeficiente de correlación
        correlacion = np.corrcoef(xi,yi)[0,1]
        
        # minimos cuadrados
        ecuacion0 = girni.linealiza_lstsq(xi,yi)

        # selecciona atipicos
        atipico_std = analiza[unabaliza]['atipico_std']
        alpha    = ecuacion0['alpha']
        beta     = ecuacion0['beta']
        fdist0   = lambda d: -10*alpha*(np.log10(d))+beta
        yi0      = fdist0(xi)
        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]

        unintervalo = 'r0' # todos
        # para exportar hacia archivo o gráfica
        ecuacion[unabaliza] = {unintervalo: ecuacion0}
        ecuacion[unabaliza][unintervalo]['correlacion'] = correlacion
        
        eq_graf[unabaliza]  = {unintervalo: {'xi_graf':xi,
                                         'yi_graf':yi,
                                         'etiqueta': par_etiqueta,
                                         'linea' : yi0,
                                         'atipicos':[xi0_e,yi0_e],
                                         'atip_etiq': etiq0_e}
                               }
        
        # Intervalos radiales en sector
        intervalo = [np.min(xi),np.max(xi)]
        frontera = analiza[unabaliza]['frontera']
        if len(frontera)>0:
            # revisar si frontera esta dentro intervalo
            frontera = np.array(frontera, dtype=float)
            revisar  = (frontera>=np.min(xi)) & (frontera<=np.max(xi))
            enintervalo = list(frontera[revisar])
            intervalo.extend(enintervalo)
            intervalo = np.array(intervalo)
            ordenar   = np.argsort(intervalo)
            intervalo = intervalo[ordenar]
        n_intervalo = len(intervalo)
        
        # analizar cada subintervalo
        p_inicio = 0
        p_desde  = 0
        p_amplia = analiza[unabaliza]['p_amplia']
        atipIntv_std = analiza[unabaliza]['atipInterv_std']
        for i in range(0,n_intervalo-1,1):
            i_eq = 'r'+str(i+1)  # indice en texto

            # puntos en subintervalo [a,b]
            a = intervalo[i]
            b = intervalo[i+1]
            subintervalo = (xi >= a) & (xi <= b)  
            xi_sub = xi[subintervalo]
            yi_sub = yi[subintervalo]
            n_sub  = len(xi_sub)
            etiq_sub = par_etiqueta[p_inicio:p_inicio + n_sub]

            # amplia sub-intervalo, mejora intersecta rectas
            detras    = p_inicio
            retrocede = detras
            if detras > p_amplia:
                retrocede = p_amplia
            delante = n_xi - (p_inicio+n_sub)# -1)
            avanza  = delante
            if delante >= p_amplia:
                avanza = p_amplia
            p_desde  = p_inicio - retrocede
            p_hasta  = p_inicio + (n_sub) + avanza
            p_inicio = p_inicio + (n_sub-1)
            
            # subintervalo, amplia puntos
            xi_a = xi[p_desde:p_hasta]
            yi_a = yi[p_desde:p_hasta]
            etiq_a = par_etiqueta[p_desde:p_hasta]
            # coeficiente de correlación
            correlacion1 = np.corrcoef(xi_a,yi_a)[0,1]

            # analiza subintervalo
            ecuacion1 = girni.linealiza_lstsq(xi_a,yi_a)
            ecuacion[unabaliza][i_eq] = ecuacion1
            ecuacion[unabaliza][i_eq]['correlacion'] = correlacion1  
            # atipicos del subintervalo extendido
            alpha  = ecuacion1['alpha']
            beta   = ecuacion1['beta']
            fdist1 = lambda d: -10*alpha*(np.log10(d))+beta
            yi1  = fdist1(xi_a)
            dyi1std = ecuacion1['error_std']
            atipico_std = analiza[unabaliza]['atipInterv_std'][i]
            dyi1 = yi_a - yi1
            atipicos = np.zeros(len(xi_a),dtype=bool)
            if np.abs(dyi1std) > casicero:
                atipicos = np.abs(dyi1) >= dyi1std*atipico_std
            xi1_e = xi_a[atipicos]
            yi1_e = yi_a[atipicos]
            etiq1_e = etiq_a[atipicos]

            # para gráfica, atipicos sin extender puntos 
            atipicos_sub = (xi1_e >= a) & (xi1_e<=b)
            xi_sub1_e = xi1_e[atipicos_sub]
            yi_sub1_e = yi1_e[atipicos_sub]
            etiq_sub1e = etiq1_e[atipicos_sub]
            eq_graf[unabaliza][i_eq] = {'atipicos': [xi_sub1_e,yi_sub1_e],
                                        'atip_etiq':etiq_sub1e}
            

            # subintervalo sin atipicos
            if len(xi1_e)>0:
                atipicoNo = np.abs(dyi1) <= dyi1std*atipico_std
                xi2 = xi_a[atipicoNo]
                yi2 = yi_a[atipicoNo]
                etiq2 = etiq_a[atipicoNo]
                # coeficiente de correlación
                correlacion2 = np.corrcoef(xi2,yi2)[0,1]
                ecuacion2 = girni.linealiza_lstsq(xi2,yi2)

                # actualiza ecuación sin atipicos intervaloy
                intervalox = ecuacion1['intervalox']
                ecuacion2['intervalox'] = intervalox.copy()
                alpha = ecuacion2['alpha']
                beta  = ecuacion2['beta']
                fdist = lambda d: -10*alpha*(np.log10(d))+beta
                intervaloy = fdist(intervalox)
                ordenar   = np.argsort(intervaloy)
                intervaloy = list(intervaloy[ordenar])
                ecuacion2['intervaloy'] = intervaloy
                ecuacion[unabaliza][i_eq] = ecuacion2
                #ecuacion[unabaliza][i_eq]['correlacion1'] = correlacion1
                ecuacion[unabaliza][i_eq]['correlacion'] = correlacion2
            
        # Revisa frontera entre subintervalos,
        # usa intersección de rectas como nueva frontera
        interv_calc = np.copy(intervalo)
        if len(intervalo) >2 and intersectar==1 : 
            for i in range(0,n_intervalo-2,1):
                ai = 'r'+str(i+1)
                bi = 'r'+str(i+2)
                ma = ecuacion[unabaliza][ai]['alpha']
                ba = ecuacion[unabaliza][ai]['beta']
                mb = ecuacion[unabaliza][bi]['alpha']
                bb = ecuacion[unabaliza][bi]['beta']

                # punto de intersección o cruce
                cruzanx = 10**((bb-ba)/(10*(mb-ma)))
                dfrontera = frontera-cruzanx
                # cruce dentro de intervalo de ecuacion
                if cruzanx > intervalo[-1]:
                    cruzanx = intervalo[-1]
                if cruzanx < intervalo[0]:
                    cruzanx = intervalo[0]
                interv_calc[i+1] = cruzanx
            
        # para grafica evalua cada subintervalo sin atipicos
        n_interv_calc = len(interv_calc)
        for i in range(0,n_interv_calc-1,1):
            i_eq = 'r'+str(i+1)
            a = interv_calc[i]
            b = interv_calc[i+1]
            subintervalo = (xi >= a) & (xi <= b)
            xi_sub = xi[subintervalo]
            yi_sub = yi[subintervalo]
            xi_graf = np.copy(xi[subintervalo])
            if not(a in xi_sub):
                xi_graf = np.concatenate(([a],xi_graf),axis=0)
            if not(b in xi_sub):
                xi_graf = np.concatenate((xi_graf,[b]),axis=0)
            
            # Evalua subintervalo con la ecuacion sin atipicos
            alpha = ecuacion[unabaliza][i_eq]['alpha']
            beta  = ecuacion[unabaliza][i_eq]['beta']
            fdist = lambda d: -10*alpha*(np.log10(d))+beta
            yi1_sub = fdist(xi_sub)
            yi_graf = fdist(xi_graf)
            eq_graf[unabaliza][i_eq]['xi_graf'] = xi_graf
            eq_graf[unabaliza][i_eq]['yi_graf'] = yi_graf
            a = np.round(np.min([xi_graf]),precision)
            b = np.round(np.max([xi_graf]),precision)
            ecuacion[unabaliza][i_eq]['intervalox'] = [a,b]
            ay = np.round(np.min([yi_graf]),precision)
            by = np.round(np.max([yi_graf]),precision)
            ecuacion[unabaliza][i_eq]['intervaloy'] = [ay,by]            
             
# SALIDA
for unabaliza in ecuacion:
    print('baliza: ',unabaliza)
    for unaecuacion in ecuacion[unabaliza]:
        unintervalo  = ecuacion[unabaliza][unaecuacion]['intervalox']
        unintervaloy = ecuacion[unabaliza][unaecuacion]['intervaloy']
        error_medio  = ecuacion[unabaliza][unaecuacion]['error_medio']
        error_std    = ecuacion[unabaliza][unaecuacion]['error_std']
        eq_latex     = ecuacion[unabaliza][unaecuacion]['eq_latex']
        errorx_medio = ecuacion[unabaliza][unaecuacion]['errorx_medio']
        errorx_std   = ecuacion[unabaliza][unaecuacion]['errorx_std']
        correlacion  = ecuacion[unabaliza][unaecuacion]['correlacion']

        print('  intervalo: ',unaecuacion)
        print('    ' + eq_latex)
        print('   ','intervalox: ',np.round(unintervalo,precision))
        print('   ','intervaloy: ',np.round(unintervaloy,precision))
        print('    correlación: ',np.round(correlacion,precision))
        print('    |error_rssi| promedio: ',np.round(error_medio,precision),
              ' , std:',np.round(error_std,precision))
        print('    |error_dist| promedio: ',np.round(errorx_medio,precision),
              ' , std:',np.round(errorx_std,precision))
    print()

# salida hacia archivo
with open(arch_ecuaciones, 'w') as outfile:
    json.dump(ecuacion, outfile) 

# GRAFICA
# 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']

if tipograf=='2D':
    for unabaliza in ecuacion:
        figura,grafica = plt.subplots()
        if escala == 'log':
            grafica.set_xscale(escala,base=escalabase)

        # todos los puntos
        unintervalo = 'r0'
        xi = eq_graf[unabaliza][unintervalo]['xi_graf']
        yi = eq_graf[unabaliza][unintervalo]['yi_graf']
        etiqueta = eq_graf[unabaliza][unintervalo]['etiqueta']
        grafica.scatter(xi,yi,marker='.')
        m = len(xi)
        for i in range(0,m,1):
            grafica.annotate(etiqueta[i],
                            (xi[i],yi[i]))
            
        # linea con todos los puntos
        fdtxt = ecuacion[unabaliza][unintervalo]['eq_latex']
        yi0 = eq_graf[unabaliza][unintervalo]['linea']
        a = np.round(np.min([xi]),2)
        b = np.round(np.max([xi]),2)
        eq_texto = fdtxt+' ; ['+str(a)+','+str(b)+']'
        grafica.plot(xi,yi0,
                     label = eq_texto,
                     linestyle='dotted')

        # lineas por cada subintervalo
        eq_interv = list(ecuacion[unabaliza].keys())
        eq_interv.pop(0)
        n_intervalo = len(eq_interv)

        for i_eq in eq_interv:
            fdtxt   = ecuacion[unabaliza][i_eq]['eq_latex']
            xi_graf = eq_graf[unabaliza][i_eq]['xi_graf']
            yi_graf = eq_graf[unabaliza][i_eq]['yi_graf']
            a = np.round(np.min([xi_graf]),precision)
            b = np.round(np.max([xi_graf]),precision)
            eq_texto = fdtxt+' ; ['+str(a)+','+str(b)+']'
            grafica.plot(xi_graf,yi_graf,
                         label = eq_texto)

            # atipicos marcados en subintervalo
            [xi1_e,yi1_e] = eq_graf[unabaliza][i_eq]['atipicos']
            etiq1_e = eq_graf[unabaliza][i_eq]['atip_etiq']
            grafica.scatter(xi1_e,yi1_e, color='red')
            m = len(etiq1_e)
            for i in range(0,m,1):
                grafica.annotate(etiq1_e[i],
                                (xi1_e[i],yi1_e[i]),
                                 color='red')
        
            # lineas de frontera
            grafica.axvline(a, color='lightblue')
            valor_frontera = str(np.round(a,precision))
            grafica.annotate(valor_frontera,
                             (a,np.max([yi,yi0])),
                             color='lightblue')
            grafica.axvline(b, color='lightblue')
            valor_frontera = str(np.round(b,precision))
            grafica.annotate(valor_frontera,
                             (b,np.max([yi,yi0])),
                             color='lightblue') 
        
        # etiquetas y títulos
        grafica.legend()
        grafica.set_ylabel(medida+'_'+modo)
        grafica.set_xlabel('distancia')
        grafica.grid(True,linestyle='dotted',
                     axis='x', which='both')
        
        untitulo = unabaliza+': '+medida+'_'+modo + ' vs distancia'
        grafica.set_title(untitulo)
        
        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

2. Captura datos. Modelo de pérdidas en propagación

El esquema prototipo para el modelo de pérdidas de señal en propagación, se registra los niveles de señal Rssi a diferentes distancias, condiciones y escenarios.

El esquema básico consta de tres balizas y un dispositivo de captura de datos. Cada baliza emite una señal a intervalos de tiempo que al ser captadas por el dispositivo permite obtener los valores Rssi y SNR. Dado que el objetivo final es revisar un modelo de localización, se requieren al menos tres balizas en el esquema.

El modelo de propagación require realizar mediciones en varios puntos a diferentes distancias que en campus presenta diferentes entornos.

Las mediciones de Rssi en un mismo punto no son constantes, por lo que se registran «muchos» valores, al menos 100 lecturas para estimar el comportamiento y usar valores representativos como el promedio y la desviación estándar..

La captura de los datos  se simplifica usando un dispositivo LoRa y un computador portátil que almacena los datos en un archivo txt.

El dispositivo LoRa se programa para tomar lecturas de Rssi y SNR de cada baliza e inmediatamente enviarlos por la conexión serial/USB, al mismo tiempo que el dispositivo se alimenta de energía por el puerto USB.

Desde un computador portátil se leen los datos capturados por el dispositivo y leídos por la conexión USB  y se almacenan en un archivo.txt. La lectura de datos serial/USB se realiza con un algoritmo en Python, que secuencialmente añade datos al archivo.txt para posterior procesamiento.

Las siguientes secciones describen lo realizado para:

Captura datos. Balizas y dispositivo – configuración algoritmo.ino

Captura datos. USB-Serial a un archivo.txt con Python

4.3 Rssi vs Distancia. Linealiza POR intervalos

En áreas extensas de medición donde existen diferentes ambientes o entornos como vegetación en una parte y edificios en otros, el resultado de básico de un solo intervalo puede mejorarse al utilizar subintervalos para cada entorno.

El cambio de entorno forma una frontera, observando la distancia al gateway (baliza) se la toma el valor como punto de partida para estimar la linealización de cada sub-intervalo. Se determinan las ecuaciones para cada intervalo, para calcular los puntos de intersección de cada linea teniendo como resultado una frontera entre diferentes ambientes calculada a partir de las mediciones realizadas.

baliza: ‘gtwFIEC’

Por ejemplo, para ‘gtwFIEC’ y su entorno mostrado en la imagen, en los puntos cercanos existe un ambiente con principalmente vegetación, que luego en puntos más alejados hay ambientes urbanos con edificios de aulas y administrativos.

Partiendo del ‘gtwFIEC’ se estima una frontera inicial a 176 metros para usar dos subintervalos y realizar la linealización.

 

Los resultados que se obtienen son:

Un intervalo

Rssi(d) = -10(4.908)log_{10}(d)+1.406 52.54 \le d \le 397.15

|error_rssi| promedio: 4.84 ,  std: 5.56

|error_dist| promedio: 41.43 ,  std: 54.09

Dos intervalos

Rssi_0(d) = -10(4.263)log_{10}(d)+(-11.269) 52.54 \le d \le 166.14

|error_rssi| promedio: 2.92 , std: 3.33

|error_dist| promedio: 20.19 ,  std: 24.67

Rssi_1(d) = -10(6.09)log_{10}(d)+(29.295) 166.14 \le d \le 397.15

|error_rssi| promedio: 2.71 , std: 3.17

|error_dist| promedio: 25.03 ,  std: 29.93

baliza: ‘gtwFCNM’

Para el caso ‘gtwFCNM’ el área de interés con vegetación se encuentra a partir de los 230 metros, tomado como valor inicial de frontera.


Un intervalo

Rssi(d) = -10(5.403)log_{10}(d)+(8.423) 27.74 \le d \le 364.71

|error_rssi| promedio: 4.59 , std: 5.48

|error_dist| promedio: 41.33 ,  std: 48.24

Dos intervalos

Rssi_0(d) = -10(4.795)log_{10}(d)+(-3.65) 27.74 \le d \le 238.81

|error_rssi| promedio: 5.96 std: 6.51

|error_dist| promedio: 52.15 ,  std: 64.38

Rssi_1(d) = -10(9.027)log_{10}(d)+(96.989) 238.81 \le d \le 364.71

|error_rssi| promedio: 2.81 , std: 3.43

|error_dist| promedio: 20.26 ,  std: 25.06

baliza: ‘gtwRECT’

Un intervalo

Rssi(d) = -10(5.054)log_{10}(d)+(11.823) 138.15 \le d \le 451.42

|error_rssi| promedio: 3.11 , std: 3.97

|error_dist| promedio: 42.93 ,  std: 51.69

Dos intervalos

Rssi_0(d) = -10(4.4443)log_{10}(d)+(-3.906) 138.15 \le d \le 358.99

|error_rssi| promedio: 1.38 , std: 1.77

|error_dist| promedio: 18.99 ,  std: 23.13

Rssi_1(d) = -10(5.905)log_{10}(d)+(33.414) 358.99 \le d \le 451.42

|error_rssi| promedio: 1.35 , std: 1.53

|error_dist| promedio: 21.17 ,  std: 24.68


Parámetros

Los parámetros para el algoritmo se establecen en el diccionario ‘analiza’. frontera es el vector donde se indica los puntos de corte del intervalo bajo estudio. Si frontera es vacio [], se asume que se analiza todo el intervalo.

Los parámetros usados en el algoritmo para obtener los resultados presentados corresponden a:

# Analizar por segmentos
analiza = {'gtwRECT':{'analizar'  : 1,
                      'atipico_std' : 1,
                      'frontera'    : [320],
                      'atipInterv_std': [1,1],
                      'p_amplia': 4,
                      'grp' : ['RECT','FIEC'],
                      'tip' : ['punto'],
                      'LOS' : [1]},
           'gtwFIEC':{'analizar'  : 1,
                      'atipico_std' : 1,
                      'frontera'    : [190], 
                      'atipInterv_std': [1,1],
                      'p_amplia': 4,
                      'grp' : ['FIEC','FCNM'],
                      'tip' : ['punto'],
                      'LOS' : [1,0]},
           'gtwFCNM':{'analizar'   : 1,
                      'atipico_std' : 1,
                      'frontera'   : [235.0],
                      'atipInterv_std': [2,2],
                      'p_amplia': 4,
                      'grp' : ['FIEC','FCNM'],
                      'tip' : ['punto'],
                      'LOS' : [1,0]}
           }

En la siguiente sección se detalla el algoritmo usado para generar las ecuaciones en cada segmento, las ecuaciones son usadas para estimar la ubicación relativa a cada baliza.

4. Rssi vs distancia. Linealización y gráfica

La relación de Rssi y distancia se puede observar usando los puntos en una gráfica 2D. El modelo teórico básico de larelación en espacio libre considera linealizar usando log10(distancia) para el eje x.

Método Empírico para ecuación Rssi(distancia)

La ecuación empírica con la que se realiza la primera estimación se modela como:

RSSI(d) = -10 \alpha \log_{10} (d) + P_{0} 1<d<L

1. Gráfica de todos los puntos en gtwFIEC

Para observar el resultado integral de todos los puntos medidos en las áreas de vegetación de FIEC y edificios FCNM realiza un primer modelo de linealización con todos los puntos. Para el eje x se usa log10(d)

Como ejemplo se presenta la gráfica para la baliza ‘d2’ etiquetada como gtwFIEC.

La linealización se realiza usando el método de los mínimos cuadrados (least square).

Al tomar como cota de error una desviación estándar por encima y por debajo de la línealización, se pueden discriminar los datos más alejados como atípicos. Se puede realizar nuevamente la linealización sin considerar los datos atípicos resultando un aumento en el valor de la pendiente α.

Al observar los datos de la gráfica que se encuentran fuera de la banda de una desviación estándar (σ) se tiene que se encuentran distribuidos en grupos: izquierda, centro y derecha. La observación lleva a considerar usar dos intervalos para realizar linealización, que es acorde los entornos de medición;  los primeros puntos con etiqueta FIEC se realizaron en zona de vegetación, y los etiquetados como FCNM se realizaron en zona de edificios de aulas y administrativos.

El análisis de los datos se divide en dos secciones:

  • Análisis con un solo entorno, un intervalo para eje distancia
  • Análisis con dos entornos, dos intervalos para eje distancia

 

 

 

 

 

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

A 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:
    
s3Eva_IT2018_T1 Intersección de dos círculos
'''
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)