Coordenadas – Rayo reflejado en perfil por tramos

Referencia: Optica – Rayo reflejado en plano inclinado

Para el análisis de rayo reflejado en un perfil de terreno, se realiza calculando el punto reflejado para cada tramo del perfil usando como base el algoritmo de rayo reflejado en plano inclinado.

Instrucciones en Python

# rayo incidente y reflejado
# en perfil por tramos
# blog.espol.edu.ec/girni
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# posición de antenas tx, rx
xi = [1,11]
yi = [4,2]

# perfil suelo
sx = [1,4,8,11]
sy = [1,0.5,0.4,1]

# PROCEDIMIENTO
def reflejotramo(xi,yi,sx,sy):
    # posición de antenas tx, rx
    xa,xb = xi
    ya,yb = yi
    # perfil suelo
    xas,xbs = sx
    sa ,sb  = sy
    # pendiente de suelo
    ms = (sb-sa)/(xbs-xas)
    bs = sa-ms*xas
    # punto de reflejo
    numerador   = bs*(xa+xb)-(xa*yb+xb*ya)+2*ms*xa*xb
    denominador = 2*bs-(ya+yb)+ms*(xa+xb)
    xc = numerador/denominador
    sc = ms*xc+bs
    reflejado = []
    if xc>=xas and xc<=xbs:
        reflejado = [xc,sc]
    return(reflejado)

# analiza para cada tramo del suelo
n = len(sx)
xc = [] ; sc = []
for k in range(0,n-1,1):
    sxk = sx[k:k+2] # un tramo de suelo
    syk = sy[k:k+2]
    reflejado = reflejotramo(xi,yi,sxk,syk)
    if len(reflejado)>0:
        xc.append(reflejado[0])
        sc.append(reflejado[1])

# SALIDA
print('puntos reflejados:')
print('xc:',xc)
print('sc',sc)

# GRAFICA
# antenas puntos en el plano
plt.scatter([xi[0],xi[1]],
            [yi[0],yi[1]],
            label='antenas')
plt.plot([xi[0],xi[1]],
         [yi[0],yi[1]],
         color='green',
         linestyle='dashed',
         label='directo')
# analiza para cada tramo del suelo
nc = len(sc)
plt.scatter(xc,sc,color='orange',
            label='punto reflejo')
# lineas de rayos
for k in range(0,nc,1):
    # lado izquierdo
    plt.plot([xi[0],xc[k]],
             [yi[0],sc[k]],
             color='blue',
             linestyle='dashdot')
    # lado derecho
    plt.plot([xc[k],xi[1]],
             [sc[k],yi[1]],
             color='blue',
             linestyle='dashdot')
    # etiquetas anotadas
    #plt.annotate(' reflejo',[xc,sc])

# perfil de suelo
plt.plot(sx,sy,color='brown',
         label='perfil suelo')

# etiquetas
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('reflexión de rayos en suelo por tramos')
plt.grid()

plt.show()

7. LoRaWan – Correlación de medidas vs estación meteorológica

La matriz de correlación entre los parámetros de propagación observados y las variables de clima de una estación meteorológica se puede procesar con Python.

Se dispone de las descripciones estadísticas de Rssi obtenidas en cada punto sobre una línea de propagación usadas para estimar el modelo de propagación para esa ruta.

Por otra parte también se dispone de las lecturas de una estación meteorológica cercana al lugar de mediciones.

La primera parte del proceso consiste en leer los archivos en diferentes tablas, para seguir un procedimiento donde se emparejan los valores de las variables a usar. El emparejamiento se realiza buscando la fecha y hora mas cercana usando la instrucción de Pandas dataframe.

    fechap = tablaPt['publishedAt'][i]
    donde  = tablaEM['fecha'].searchsorted(fechap)

La matriz de correlación  se calcula usando la librería Pandas para mantener un formato de columnas de lectura sencilla. Los datos en columnas a usar se agrupan en una tabla «mediciones» para realizar la llamada a corr().

correlacion_matriz = mediciones.corr()

la tabla de correlación que se obtiene:

Matriz de correlación: m_MD10
             rssi_up  rssi_down      TEMP  Humidity  Solar_rad  Bar_press.  Rainfall
rssi_up     1.000000   0.062009 -0.103856  0.050208  -0.079344    0.172215       NaN
rssi_down   0.062009   1.000000 -0.058020 -0.004584   0.002753    0.158011       NaN
TEMP       -0.103856  -0.058020  1.000000 -0.972939   0.841263   -0.673707       NaN
Humidity    0.050208  -0.004584 -0.972939  1.000000  -0.819612    0.618423       NaN
Solar_rad  -0.079344   0.002753  0.841263 -0.819612   1.000000   -0.361084       NaN
Bar_press.  0.172215   0.158011 -0.673707  0.618423  -0.361084    1.000000       NaN
Rainfall         NaN        NaN       NaN       NaN        NaN         NaN       NaN

La matriz de correlación se aplica para cada punto.

Instrucciones en Python

# Archivos de sensores y estación meteorologica
# correlacion entre medida_Punto y sensor_EM
import numpy as np
import pandas as pd
import datetime as dt
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter, DayLocator

# INGRESO
# tabla punto
archivoPunto = "data_m_MD10.csv"
medida_Punto = ["rssi_up","rssi_down"]
medida_Punto_u = ["(dBm)","(dBm)"]

# tabla estacion meteorologica
archivoEM = "2021_10a12EstMeteorologica.csv"
sensor_EM = ["TEMP","Humidity","Solar_rad",
             "Bar_press.","Rainfall"]
sensor_EM_u = ["(°C)","(%)","(kW)",
               "(hPa)","(mm)"]

# carpeta de resultados
carpeta_rsm = "correlacion_rsm"

# PROCEDIMIENTO
# tabla punto
tablaPt = pd.read_csv(archivoPunto)
tablaPt = tablaPt.drop(columns='Unnamed: 0')
tablaPt = pd.DataFrame(tablaPt)
tablaPt_n = len(tablaPt)
fechaformatoPt = "%Y-%m-%d %H:%M:%S.%f"
tablaPt['publishedAt'] = pd.to_datetime(tablaPt['publishedAt'],
                                      format=fechaformatoPt)
# tabla estacion meteorologica
tablaEM = pd.read_csv(archivoEM, sep=';',decimal=',')
tablaEM = pd.DataFrame(tablaEM)
tablaEM_n = len(tablaEM)
fechaformatoEM = "%d/%m/%Y %H:%M:%S"
tablaEM['fecha'] = tablaEM['Date']+' '+tablaEM['Time']
tablaEM['fecha'] = pd.to_datetime(tablaEM['fecha'],
                                  format=fechaformatoEM)

# intersecta tablas: punto con estacion meteorologica
# columnas vacias
for j in range(0,len(sensor_EM),1):
    tablaPt[sensor_EM[j]] = np.NAN
# busca fecha y hora mas cercana
for i in range(0,tablaPt_n-1,1):
    fechap = tablaPt['publishedAt'][i]
    donde  = tablaEM['fecha'].searchsorted(fechap)
    if donde<tablaEM_n:
        for j in range(0,len(sensor_EM),1):
            sensor_EMvalor = tablaEM[sensor_EM[j]][donde]
            tablaPt.at[i,sensor_EM[j]] = sensor_EMvalor

# Matriz de Correlación
columnas = medida_Punto.copy()
columnas.extend(sensor_EM)
mediciones = tablaPt[columnas].copy()
correlacion_matriz = mediciones.corr()

# Para graficas
ti  = tablaPt['publishedAt']
medida_graf = {}
untitulo = archivoPunto.split('.')[0][5:]
for i in range(0,len(medida_Punto),1):
    xi = tablaPt[medida_Punto[i]]
    xiavg = np.around(np.average(xi),1)
    xistd = np.around(np.std(xi),2)
    etiq = medida_Punto[i]+' '+str(xiavg)
    etiq = etiq+' +/- '+str(xistd)
    medida_graf[i] = {'medida':xi,'promedio':xiavg,
                      'std':xistd,'etiqueta':etiq}
# SALIDA
print('Matriz de correlación:', untitulo)
print(correlacion_matriz)

# Guarda archivo de matriz de correlación
unresumen = carpeta_rsm+'/'+untitulo+'_correlacion.csv'
tablacorr = correlacion_matriz.round(decimals=4)
tablacorr.to_csv(unresumen)

# Grafica  variable vs tiempo --------
colores = ['blue','orange','green','magenta']
for j  in range(0,len(sensor_EM),1): # cada sensor_EM
    graf_ylim =[]
    fig_Pt, graf_EM = plt.subplots(2,1)
    for i in range(0,len(medida_Punto),1): # cada medida_Punto
        # medida del punto
        graf_EM[i].scatter(ti,medida_graf[i]['medida'],
                           marker = '.',color=colores[i],
                           label=medida_graf[i]['etiqueta'])
        # estacion meteorologica
        eje12 = graf_EM[i].twinx()  # eje de izquierda
        yi = tablaPt[sensor_EM[j]]
        eje12.scatter(ti,yi, marker = '.',
                      color='lightgreen',label=sensor_EM[j])
        # etiquetas
        graf_EM[i].xaxis.set_major_formatter(DateFormatter('%H:%M'))
        etiq_izq = medida_Punto[i]+' '+medida_Punto_u[i]
        graf_EM[i].set_ylabel(etiq_izq, color=colores[i])
        etiq_der = sensor_EM[j]+' '+sensor_EM_u[j]
        eje12.set_ylabel(etiq_der, color='green')
        graf_EM[i].legend()
        graf_ylim.extend(list(graf_EM[i].get_ylim()))

    # ejey subgraficas, limites izquierdo iguales
    y_a = np.min(graf_ylim)
    y_b = np.max(graf_ylim)
    for i in range(0,len(medida_Punto),1):
        graf_EM[i].set_ylim(y_a,y_b)
    # titulos
    titulo_graf = untitulo +': '+medida_Punto[0].split('_')[0]
    titulo_graf = titulo_graf +' vs '+sensor_EM[j]
    graf_EM[0].set_title(titulo_graf)
    fig_Pt.tight_layout()
    # archivo de grafico
    arch_graf = untitulo +'_'+medida_Punto[0].split('_')[0]
    arch_graf = arch_graf +'_'+sensor_EM[j]
    plt.savefig(carpeta_rsm+'/'+arch_graf+'.png')
# plt.show()

# Grafica entre variables --------
for j  in range(0,len(sensor_EM),1): # cada sensor_EM
    fig_Pt, graf_EM = plt.subplots(2,1)
    for i in range(0,len(medida_Punto),1): # cada medida_Punto
        graf_EM[i].scatter(medida_graf[i]['medida'],
                           tablaPt[sensor_EM[j]],
                           marker = '.',color=colores[i])
        # etiquetas
        etiq_x = medida_Punto[i]+' '+medida_Punto_u[i]
        graf_EM[i].set_xlabel(etiq_x, color=colores[i])
        etiq_y = sensor_EM[j]+' '+sensor_EM_u[j]
        graf_EM[i].set_ylabel(etiq_y, color='green')
    # titulos
    titulo_graf = untitulo +': '+medida_Punto[0].split('_')[0]
    titulo_graf = titulo_graf +' vs '+sensor_EM[j]
    graf_EM[0].set_title(titulo_graf)
    fig_Pt.tight_layout()
    # archivo de grafico
    arch_graf = untitulo +'_'+medida_Punto[0].split('_')[0]
    arch_graf = arch_graf +'_'+sensor_EM[j]
    plt.savefig(carpeta_rsm+'/'+arch_graf+'.png')
plt.show()

Referencia: Correlación y lectura de archivos con Python

Correlación con Python – Ejercicio Estación Meteorológica. ESTG1003 – Procesos Estocásticos

Archivos.csv con Python – Ejercicio con gráfica de temperatura y Humedad. CCPG1001 – Fundamentos de programación

6. LoRaWan – Linealiza por segmentos

Para los diferentes entornos donde se propaga la señal, se diferencian por segmentos (fronteras) delimitados en metros desde el gateway.

El desnivel del terreno que existe aproximandamente a 90 mts desde el gateway tiene una sección de «sombra» en donde no se usarán los valores para generar el modelo. por lo que el segundo segmento no se usará para el cálculo, se indica en la variable como [1,0,1]

frontera = [0,90,125,250] # metros 
frontera_usar = [1,0,1] # usar

Usando estos parámetros y reordenando el algoritmo anterior para ser usado en cada segmento, se obtiene el siguiente resultado

Para la sección de sombra, para mantener la continuidad en la gráfica, se une el punto final del segmento anterior, y el primer punto del segmento posterior.

Los detalles de cada ecuación por segmentos son:

 ---- segmento: 0
                                            ecuacion_up                             ecuacion_down
alpha                                          2.579875                                  2.366443
beta                                         -25.528314                                -27.974302
eq_latex          $ rssi = -10(2.58)log_{10}(d)-25.53 $     $ rssi = -10(2.37)log_{10}(d)-27.97 $
intervalox                                 [1.0, 86.05]                              [1.0, 86.05]
error_medio                                     6.17513                                  6.495915
error_std                                      7.697343                                  7.230741
eqg_latex       $ d = 10^{(-25.53 - rssi)/(10(2.58))} $   $ d = 10^{(-27.97 - rssi)/(10(2.37))} $
intervaloy    [-90.48339483394834, -25.745454545454542]  [-82.15682656826569, -26.98787878787879]
errorx_medio                                  43.575249                                 35.657177
errorx_std                                    89.139552                                 50.786448

 ---- segmento: 1
                                           ecuacion_up                             ecuacion_down
alpha                                         7.068273                                  6.777292
beta                                         61.310998                                 57.364631
eq_latex         $ rssi = -10(7.07)log_{10}(d)+61.31 $     $ rssi = -10(6.78)log_{10}(d)+57.36 $
intervalox                             [86.05, 124.08]                           [86.05, 124.08]
error_medio                                        0.0                                       0.0
error_std                                          0.0                                       0.0
eqg_latex       $ d = 10^{(61.31 - rssi)/(10(7.07))} $    $ d = 10^{(57.36 - rssi)/(10(6.78))} $
intervaloy    [-86.67755403724453, -75.44247020329829]  [-84.53164627778156, -73.75907940898257]
errorx_medio                                       0.0                                       0.0
errorx_std                                         0.0                                       0.0

 ---- segmento: 2
                                            ecuacion_up                             ecuacion_down
alpha                                          4.082007                                  4.159709
beta                                          -1.212499                                  2.560248
eq_latex           $ rssi = -10(4.08)log_{10}(d)-1.21 $      $ rssi = -10(4.16)log_{10}(d)+2.56 $
intervalox                             [124.08, 238.66]                          [126.91, 235.77]
error_medio                                    3.623762                                  3.887403
error_std                                      4.452527                                  4.535259
eqg_latex        $ d = 10^{(-1.21 - rssi)/(10(4.08))} $     $ d = 10^{(2.56 - rssi)/(10(4.16))} $
intervaloy    [-105.64163822525596, -81.32270916334662]  [-99.64505119453923, -80.30278884462152]
errorx_medio                                  39.417734                                 39.102089
errorx_std                                     52.33983                                 45.251062

Instrucciones en Python

# ecuación por mínimos cuadrados de una ruta por segmentos
# usando las medidas y distancia al gateway de cada punto
# Revision 2022/07/23
# http://blog.espol.edu.ec/girni/

import numpy as np
import pandas as pd
import os
import json
import matplotlib.pyplot as plt
import girni_lorawan_lib as girni

# INGRESO
carpeta_rsm = "resultado_c_Maiz_todo"
frontera = [0,90,125,250] # metros
segmento_usar =  [1,0,1] # usar

arch_var_general = 'variables_generales.txt'

# PROCEDIMIENTO
# carga variables desde archivos
with open(arch_var_general) as linea:
    texto = linea.read()
var_gen = json.loads(texto)
globals().update(var_gen)

partes = carpeta_rsm.strip('/').split('_')
arch_nombre = partes[1]+'_'+partes[2]

def linealiza_preparadatos(carpeta_rsm,arch_coord,var_gen):
    ''' desarrolla la ecuación para un gradiente
    '''
    var_gen['movAvg_cual'] = [2,4] #cada cuantas muestras
    partes = carpeta_rsm.strip('/').split('_')
    arch_nombre = partes[1]+'_'+partes[2]

    carp_coord = var_gen['carp_coord']
    medida = var_gen['medida']
    precision = var_gen['precision']
    
    # lee coordenadas y su distancia al gateway
    dist_Gtw = girni.coordenadas_leer(arch_coord,carp_coord)

    # lista.txt de archivos a usar
    arch_lista = medida+"_"+arch_nombre+"_lista.txt"
    archivo_ruta  = carpeta_rsm + '/' + arch_lista
    archivoexiste = os.path.exists(archivo_ruta)
    if not(archivoexiste):
        puntoUsar_Colum = ['punto', 'up', 'down'] 
        puntoUsar_Tipos = {'punto' : 'object',
                           'up'   : 'int64',
                           'down' : 'int64'}
        puntoUsar = pd.DataFrame(columns = puntoUsar_Colum)
        puntoUsar = puntoUsar.astype(dtype = puntoUsar_Tipos)
        for unarchivo in os.listdir(carpeta_rsm):
            verifica = unarchivo.startswith('describe_')
            verifica = verifica and unarchivo.endswith('.csv')
            if verifica:
                partes  = unarchivo.strip('.csv').split('_')
                unpunto = partes[2]
                puntodato = {'punto' : unpunto,
                             'up' : 1,
                             'down' : 1}
                puntodato = pd.DataFrame([puntodato])
                puntoUsar = pd.concat([puntoUsar,
                                       puntodato],
                                      ignore_index = True)
        puntoUsar = puntoUsar.set_index('punto')
        puntoUsar.to_csv(archivo_ruta)
    # usa lista.txt de archivos seleccionados 1 ó 0 para enlace up,down
    if (archivoexiste):
        puntoUsar = pd.read_csv(archivo_ruta)
        puntoUsar = puntoUsar.set_index('punto')
      
    # Datos para gráfica desde lista.txt
    punto_columnas = ['codigopunto','dist_Gtw',
                      medida+'_up',medida+'_up'+'_std',
                      'usar_up',medida+'_down',
                      medida+'_down'+'_std',
                      'usar_down','dispositivo']
    punto_tipos = {'codigopunto': 'object',
                   'dist_Gtw'   : 'float64',
                   medida+'_up' : 'float64',
                   medida+'_up'+'_std': 'float64',
                   'usar_up' : 'int64',
                   medida+'_down': 'float64',
                   medida+'_down'+'_std' :'float64',
                   'usar_down' : 'int64',
                   'dispositivo': 'object'}
    punto_graf = pd.DataFrame(columns=punto_columnas)
    punto_graf = punto_graf.astype(dtype=punto_tipos)
    puntoSinDist = []
    for unarchivo in os.listdir(carpeta_rsm):
        if unarchivo.startswith('describe_'):
            codigopunto = unarchivo.strip('.csv').split('_')[2]
            
            # lectura del archivo
            unresumen  = carpeta_rsm+'/'+unarchivo
            descriptor = pd.read_csv(unresumen,index_col='Unnamed: 0')

            if (codigopunto in dist_Gtw.index):
                undato = {'codigopunto': codigopunto,
                          'dist_Gtw'   : dist_Gtw[codigopunto],
                          medida+'_up' : descriptor['rssi_up']['mean'],
                          medida+'_up'+'_std': descriptor['rssi_up']['std'],
                          'usar_up' : puntoUsar['up'][codigopunto],
                          medida+'_down': descriptor['rssi_down']['mean'],
                          medida+'_down'+'_std' :descriptor['rssi_down']['std'],
                          'usar_down' : puntoUsar['down'][codigopunto],
                          'dispositivo': descriptor['dispositivo'][0]
                          }
                undato = pd.DataFrame([undato])
                punto_graf = pd.concat([punto_graf,undato],
                                       ignore_index=True)
            else:
                puntoSinDist.append(codigopunto)
            
    punto_graf = punto_graf.set_index('codigopunto')
    punto_graf = punto_graf.sort_values('dist_Gtw' )

    return(punto_graf)


punto_graf = linealiza_preparadatos(carpeta_rsm,arch_coord,var_gen)

# Segmentos, intervalos de indices en tabla ----------
extiende   = 15 # metros
intervalos = []
intervalos_ext = []

n_sector = len(frontera)
indices = np.arange(0,len(punto_graf['dist_Gtw']),1)
for i in range(0,n_sector-1,1):
    desde = frontera[i]
    hasta = frontera[i+1]
    
    # intervalos, indices [a,b]
    sector = (punto_graf['dist_Gtw']>=desde)
    sector = sector & (punto_graf['dist_Gtw']<hasta)
    unintervalo = indices[sector]
    a = np.min(unintervalo)
    if i>0:
        a = a - 1
        sector[a] = True
    b = np.max(unintervalo)
    intervalos.append([a,b])
    punto_graf['sector_'+str(i)] = sector
    
    # intervalos extendidos
    if (desde-extiende)>0:
        desde = desde - extiende
    if (hasta + extiende) < np.max(punto_graf['dist_Gtw']):
        hasta = hasta + extiende
    sector_ext = (punto_graf['dist_Gtw']>=desde)
    sector_ext = sector_ext & (punto_graf['dist_Gtw']<=hasta)
    unintervalo = indices[sector_ext]
    a = np.min(unintervalo)
    b = np.max(unintervalo)
    intervalos_ext.append([a,b])
    

def linealiza_segmento(punto_graf,var_gen,segmento=-1):
    '''estima la ecuacion para cada segmento,
       usada solo cuando segmento >=0 [0,1,2,... ]
    '''
    precision = var_gen['precision']
    resultado  = {}
    xi = np.array(punto_graf['dist_Gtw'])
    # enlace up/down
    enlaces =['up','down']
    for unenlace in enlaces:
        usar = punto_graf['usar_'+unenlace]
        if segmento>=0:
            sector = punto_graf['sector_'+str(segmento)]
            usar   = usar*sector
        yi = np.array(punto_graf[medida+'_'+unenlace])
        # ecuacion
        eq = girni.linealiza_lstsq(xi[usar==1],
                                   yi[usar==1],
                                   digitos = precision)
        alpha = eq['alpha']
        beta  = eq['beta']
        fd = lambda d: -10*alpha*(np.log10(d))+beta
        # puntos de ecuacion en tabla
        punto_graf['fi_'+unenlace] = fd(xi)
        resultado['ecuacion_'+unenlace] = eq
    resultado['tabla'] = punto_graf
    return (resultado)

# segmentos a usar primero y último  ----------
primero   = segmento_usar.index(1)
invertido = segmento_usar.copy()
invertido.reverse()
ultimo = len(segmento_usar) - invertido.index(1) - 1

# ecuaciones para todos los segmentos
eq_intervalo = {}
tabla = {}
for i in range(0,n_sector-1,1):
    resultado = linealiza_segmento(punto_graf,var_gen,i)
    ecuacion = resultado.copy()
    ecuacion.pop('tabla')
    ecuacion = pd.DataFrame(ecuacion)
    eq_intervalo[i] = ecuacion
    tabla[i] = resultado['tabla'].copy()
        
# intervalos de sombra corrige anteriores
for i in range(0,n_sector-1,1):
    condicion = i>primero and i<ultimo
    if not(segmento_usar[i]) and i>0 and condicion:
        sector = punto_graf['sector_'+str(i)]
        # extremos [a,b]
        a = intervalos[i-1][1]
        b = intervalos[i][1]
        # extremo izquierdo a
        xa_up = tabla[i]['dist_Gtw'][a]
        ya_up = tabla[i-1]['fi_up'][a]
        ya_down = tabla[i-1]['fi_down'][a]
        # extremo derecho b
        xb_up = tabla[i]['dist_Gtw'][b]
        yb_up = tabla[i+1]['fi_up'][b]
        yb_down = tabla[i+1]['fi_down'][b]
        eq_up = girni.linealiza_lstsq([xa_up,xb_up],
                                      [ya_up,yb_up],
                                      digitos = precision)
        eq_down = girni.linealiza_lstsq([xa_up,xb_up],
                                        [ya_down,yb_down],
                                        digitos = precision)
        ecuacion = {'ecuacion_up'   : eq_up,
                    'ecuacion_down' : eq_down}
        ecuacion = pd.DataFrame(ecuacion)
        eq_intervalo[i] = ecuacion
        
# SALIDA
for i in range(0,n_sector-1,1):
    print('\n ---- segmento:',i)
    print(eq_intervalo[i])

def graf_gradiente_segmento(arch_nombre,punto_graf,ecuacion,var_gen):
    ''' grafica de gradiente
    '''
    medida = var_gen['medida']
    precision = var_gen['precision']
    escalabase = 10    # 10
    escala = 'log'

    dist = punto_graf['dist_Gtw']
    # enlace up/down
    enlaces = ['up','down']
    colorenlace = ['blue','orange']
    i=0
    for unenlace in enlaces:
        # enlace_up/down
        media = punto_graf[medida+'_'+unenlace]
        std   = punto_graf[medida+'_'+unenlace+'_std']
        media_techo = media + std
        media_piso  = media - std
        yi = np.array(media)
        fi = punto_graf['fi'+'_'+unenlace]

        usar = punto_graf['usar_'+unenlace]
        eq = ecuacion['ecuacion_'+unenlace]
        dyi_std = eq['error_std']

        for j in range(0,len(dist),1):
            uncolor = colorenlace[i]
            etiqueta = dist.keys()[j]
            if not(usar[j]):
                uncolor = 'red'
            # medida up/down +/- std
            graf[i].scatter(dist[j],media[j],
                            color=uncolor,marker='o')
            graf[i].scatter(dist[j],media_techo[j],
                            color=uncolor,marker='+')
            graf[i].scatter(dist[j],media_piso[j],
                            color=uncolor,marker='+')
            # etiquetas por punto
            graf[i].annotate(etiqueta,(dist[j],yi[j]),rotation=45)
        
        # linealizado up/down
        etiq_gradnt = eq['eq_latex']
        etiq_gradnt = etiq_gradnt+' ; ['+str(int(np.min(dist)))
        etiq_gradnt = etiq_gradnt+','+str(int(np.max(dist)))+']'
        graf[i].plot(dist,fi,label = etiq_gradnt,color=colorenlace[i])
        graf[i].plot(dist,fi+dyi_std,color=colorenlace[i],linestyle='dotted')
        graf[i].plot(dist,fi-dyi_std,color=colorenlace[i],linestyle='dotted')

        # etiquetado de grafico
        graf[i].set_ylabel(medida+'_'+enlaces[i]+' (dBm)',color=colorenlace[i])
        graf[i].legend()
        graf[i].grid(True,linestyle='dotted',axis='x', which='both')
        if i==0:
            graf[0].set_title(arch_nombre+' '+medida)
        if i==1:
            graf[1].set_xlabel('distancia')
        i = i + 1
    return(fig_gradnt)

# grafica gradientes por segmentos
# para limites en eje y
techo_up = pd.DataFrame()
piso_up  = pd.DataFrame()
techo_down = pd.DataFrame()
piso_down  = pd.DataFrame()

graf=[0,0]
fig_gradnt,graf = plt.subplots(2,1)
escalabase = 10    # 10
escala = 'log'
if escala == 'log':
    graf[0].set_xscale(escala,base=escalabase)
    graf[1].set_xscale(escala,base=escalabase)

# grafica cada segmento
for i in range(0,n_sector-1,1):
    sector = tabla[i]['sector_'+str(i)]
    punto_grupo = tabla[i][sector]

    #enlace_up
    media_up = punto_grupo[medida+'_up']
    std_up   = eq_intervalo[i]['ecuacion_up']['error_std']
    media_up_techo = media_up + std_up
    media_up_piso  = media_up - std_up
    # enlace_down
    media_down = punto_grupo[medida+'_down']
    std_down   = eq_intervalo[i]['ecuacion_up']['error_std']
    media_down_techo = media_down + std_down
    media_down_piso  = media_down - std_down
    
    techo_up = pd.concat([techo_up,media_up_techo])
    piso_up  = pd.concat([piso_up,media_up_piso])
    techo_down = pd.concat([techo_down,media_down_techo])
    piso_down  = pd.concat([piso_down,media_down_piso])

    if segmento_usar[i]:
        ecuacion = eq_intervalo[i]
        graf_gradiente_segmento(arch_nombre,punto_grupo,ecuacion,var_gen)
        
    condicion = i>primero and i<ultimo
    if not(segmento_usar[i]) and i>0 and condicion:
        a = tabla[i]['dist_Gtw'][intervalos[i][0]]
        b = tabla[i]['dist_Gtw'][intervalos[i][1]]
        
        graf[0].axvspan(a,b,alpha=0.5,color='gray')
        graf[1].axvspan(a,b,alpha=0.5,color='gray')
        
        # linea intermedia
        xi = punto_grupo['dist_Gtw']
        ecuacion = eq_intervalo[i]
        fi_up = punto_grupo['fi_up']
        fi_down = punto_grupo['fi_down']

        graf[0].scatter(xi,media_up,color='blue',marker='o')
        graf[0].scatter(xi,media_up_techo,color='blue',marker='+')
        graf[0].scatter(xi,media_up_piso,color='blue',marker='+')
        
        graf[1].scatter(xi,media_down,color='orange',marker='o')
        graf[1].scatter(xi,media_down_techo,color='orange',marker='+')
        graf[1].scatter(xi,media_down_piso,color='orange',marker='+')

        alpha_up = eq_up['alpha']
        beta_up  = eq_up['beta']
        fd_up = lambda d: -10*alpha_up*(np.log10(d))+beta_up
        fi_up = fd_up(xi)
        
        alpha_down = eq_down['alpha']
        beta_down  = eq_down['beta']
        fd_down = lambda d: -10*alpha_down*(np.log10(d))+beta_down
        fi_down = fd_down(xi)
        
        etiq = eq_up['eq_latex']
        etiq = etiq+' ; ['+str(int(a))+','+str(int(b))+']'
        graf[0].plot(xi,fi_up, label = etiq,
                     color='blue', linestyle='dashed')

        etiq = eq_down['eq_latex']
        etiq = etiq+' ; ['+str(int(a))+','+str(int(b))+']'
        graf[1].plot(xi,fi_down, label = etiq,
                     color='brown',linestyle='dashed')

# ajuste de eje y
y_min = np.min([piso_up.min(),piso_down.min()])
y_max = np.max([techo_up.max(),techo_down.max()])
if y_min<-135:
    y_min = -135
#graf[0].set_ylim(y_min,y_max)
#graf[1].set_ylim(y_min,y_max)
plt.tight_layout()

unarchivo = carpeta_rsm+'/ecuacion_segmentos_'+arch_nombre+'.png'
fig_gradnt.savefig(unarchivo)

plt.show()

5. LoRaWan – Linealiza un intervalo

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

Las descripciones estadísticas de Rssi obtenidas en cada punto sobre una línea de propagación se usan para estimar el modelo de propagación para esa ruta.

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

La ruta de ejemplo sigue la dirección desde el gateway (Gw03) hacia una parcela de plantación de maiz (MA##) en la parte superior izquierda de la imagen. se añade la primera sección (M0##)  para cubrir todo el segmento de propagación desde el gateway.

La imagen presentada representa la ruta c_MA, la letra «c» identifica el modelo de dispositivo usado (c: cápsula, m: módulo)

Los valores de Rssi_media en cada punto, tanto los enlaces de subida (up) y bajada (down), se linealizan por el método de los mínimos cuadrados usando la distancia como variable independiente.

En cada punto Rssi se añaden las marcas ‘+’ que representan media+std y media-std.

Mejorando el modelo

Las perdidas de señal desde el gateway aumentan con la distancia, siendo más grandes al final del segmento para el canal de bajada (down), mostrando una mayor desviación estándar (std) que afectan al modelo. Para discriminar los puntos que no representan una buena cobertura, se utiliza una gráfica de desviaciones estándar, donde los valores fuera de rango no se usan en el modelo.

La selección de los puntos que incluyen en el modelo se realiza en el archivo «rssi_c_MA_lista.txt»

punto,up,down
M001,1,1
M002,1,1
M003,1,1
M004,1,1
M005,1,1
M006,1,1
M007,1,1
MA03,1,0
MA04,1,1
MA05,1,1
MA06,1,1
MA07,1,1
MA08,1,0
MA09,1,0
MA10,1,0
Ref01,1,1

En el archivo se muestra que los puntos MA03, MA08, MA09, MA10 no se usan para el enlace de bajada (down).

Se guarda el archivo de selección con la configuración, y se vuelve a ejecutar el algoritmo, y se obtiene el modelo mejorado de pérdidas de propagación:

 ecuaciones: 
                                        ecuacion_up                             ecuacion_down
alpha                                         3.01                                      2.58
beta                                        -21.37                                    -25.98
coef_correlacion                             -0.94                                     -0.93
coef_determinacion                            0.88                                      0.86
eq_latex     $ rssi = -10(3.01)log_{10}(d)-21.37 $     $ rssi = -10(2.58)log_{10}(d)-25.98 $
intervalox                           [1.0, 233.86]                             [1.0, 198.31]
error_medio                                   5.82                                      5.64
error_std                                     6.63                                      6.46
eqg_latex  $ d = 10^{(-21.37 - rssi)/(10(3.01))} $   $ d = 10^{(-25.98 - rssi)/(10(2.58))} $
intervaloy                        [-97.87, -25.75]                          [-87.68, -26.99]
errorx_medio                                 43.32                                     37.68
errorx_std                                   56.01                                     47.11

Una forma complementaria para observar los puntos «aberrantes» o fuera de rango es observar en cada punto la diferencia entre las medias del enlace de subida y el de bajada, donde se destaca que los valores se encuentran muy separados.

Las gráficas se realizan en bloques como funciones (def), a fin de activarlas o no durante el análisis detallado, además de incorporarlas a la librería girni_lorawan_lib.py

Las rutas posibles para este ejercicio se dan por medio de «carpeta_rsm», que es el directorio donde se encuentran los archivos con los detalles de cada punto.


Instrucciones en Python

Algunos parámetros se repiten entre algoritmos, por lo que se los ha separado en un archivo de texto ‘variables_generales.txt’.

El contenido de las ‘variables_generales.txt’ se muestra a continuación:

{"carpeta_DB"  : "BaseDatosJSON",
 "medida"         : "rssi",
 "medida_unidad"  : "dBm",
 "medida_normal"  : [-250,-1],
 "medida_grafica" : [-100,-60],
 "movAvg_cual"    : [8,16],
 "movAvg_color"   : ["lightgreen","orange"],
 "precision"  : 2,
 "gatewayDB"  : {"uCfr//5zvhk=":"Gw01",
                 "uCfr//4dylc=":"Gw03"},
 "arch_coord" : "coord_puntos.csv",
 "carp_coord" : "coordenadas",
 "zonaNum"    : 17,
 "zonaLetra"  : "M",
 "alturaGw"   : 8,
 "alturapunto": 1,
 "nFres"      : 1,
 "freq"       : 915,
 "muestras"   : 41,
 "guarda"     : 1
}

Con lo que las instrucciones en Python quedan como:

# Obtiene ecuación por mínimos cuadrados de una ruta
# usando las medidas y distancia al gateway de cada punto
# http://blog.espol.edu.ec/girni/

import numpy as np
import pandas as pd
import os
import json
import matplotlib.pyplot as plt
import girni_lorawan_lib as girni

# INGRESO
carpeta_rsm = "resultado_c_Maiz_todo"
arch_coord  = "coord_puntos.csv"
arch_var_general = "variables_generales.txt"

# PROCEDIMIENTO
# carga variables desde archivos
with open(arch_var_general) as linea:
    texto = linea.read()
var_gen = json.loads(texto)
globals().update(var_gen)

var_gen['movAvg_cual'] = [2,4] #cada cuantas muestras

partes = carpeta_rsm.strip('/').split('_')
arch_nombre = partes[1]+'_'+partes[2]

def coordenadas_leer(arch_coord,carp_coord):
    ''' lista las coordenadas desde un archivo
        ubicado en la carpeta
    '''
    if len(carp_coord)>0:
        carp_coord = carp_coord+'/'
    archivo_ruta = carp_coord+arch_coord
    arch_coord_existe = os.path.exists(archivo_ruta)

    if arch_coord_existe:
        puntos = pd.read_csv(archivo_ruta,index_col='punto')
        # busca gateway y ubicación de punto observado
        Gw_nombre = ''
        for unpunto in puntos.index:
            if unpunto.startswith('Gw'):
                Gw_nombre = unpunto
        dist_Gtw = puntos['dist_'+Gw_nombre]   
    else:
        print(' ERROR: NO se encuentra el archivo...')
        print('        revise el ruta/nombre de archivo. ')
    return(dist_Gtw)

dist_Gtw = coordenadas_leer(arch_coord,carp_coord)

# crea lista archivos a usar si no existe
arch_lista = medida+"_"+arch_nombre+"_lista.txt"
archivo_ruta  = carpeta_rsm + '/' + arch_lista
archivoexiste = os.path.exists(archivo_ruta)
if not(archivoexiste):
    puntoUsar_Colum = ['punto', 'up', 'down'] 
    puntoUsar_Tipos = {'punto' : 'object',
                       'up'   : 'int64',
                       'down' : 'int64'}
    puntoUsar = pd.DataFrame(columns = puntoUsar_Colum)
    puntoUsar = puntoUsar.astype(dtype = puntoUsar_Tipos)
    for unarchivo in os.listdir(carpeta_rsm):
        verifica = unarchivo.startswith('describe_')
        verifica = verifica and unarchivo.endswith('.csv')
        if verifica:
            partes  = unarchivo.strip('.csv').split('_')
            unpunto = partes[2]
            puntodato = {'punto' : unpunto,
                         'up' : 1,
                         'down' : 1}
            puntodato = pd.DataFrame([puntodato])
            puntoUsar = pd.concat([puntoUsar,
                                   puntodato],
                                  ignore_index = True)
    puntoUsar = puntoUsar.set_index('punto')
    puntoUsar.to_csv(archivo_ruta)

if (archivoexiste):
    puntoUsar = pd.read_csv(archivo_ruta)
    puntoUsar = puntoUsar.set_index('punto')
  
# Analiza cada punto en lista
punto_columnas = ['codigopunto','dist_Gtw',
                  medida+'_up',medida+'_up'+'_std',
                  'usar_up',medida+'_down',
                  medida+'_down'+'_std',
                  'usar_down','dispositivo']
punto_tipos = {'codigopunto': 'object',
               'dist_Gtw'   : 'float64',
               medida+'_up' : 'float64',
               medida+'_up'+'_std': 'float64',
               'usar_up' : 'int64',
               medida+'_down': 'float64',
               medida+'_down'+'_std' :'float64',
               'usar_down' : 'int64',
               'dispositivo': 'object'}
punto_graf = pd.DataFrame(columns=punto_columnas)
punto_graf = punto_graf.astype(dtype=punto_tipos)
puntoSinDist =[]
for unarchivo in os.listdir(carpeta_rsm):
    if unarchivo.startswith('describe_'):
        codigopunto = unarchivo.strip('.csv').split('_')[2]
        
        # lectura del archivo
        unresumen  = carpeta_rsm+'/'+unarchivo
        descriptor = pd.read_csv(unresumen,index_col='Unnamed: 0')

        # Para Grafica
        condicion = (codigopunto in dist_Gtw.index)
        if condicion:
            undato = {'codigopunto': codigopunto,
                      'dist_Gtw'   : dist_Gtw[codigopunto],
                      medida+'_up' : descriptor['rssi_up']['mean'],
                      medida+'_up'+'_std': descriptor['rssi_up']['std'],
                      'usar_up' : puntoUsar['up'][codigopunto],
                      medida+'_down': descriptor['rssi_down']['mean'],
                      medida+'_down'+'_std' :descriptor['rssi_down']['std'],
                      'usar_down' : puntoUsar['down'][codigopunto],
                      'dispositivo': descriptor['dispositivo'][0]
                      }
            undato = pd.DataFrame([undato])
            punto_graf = pd.concat([punto_graf,undato],
                                   ignore_index=True)
        else:
            puntoSinDist.append(codigopunto)
punto_graf = punto_graf.set_index('codigopunto')
punto_graf = punto_graf.sort_values('dist_Gtw' )

# selecciona datos
dist = punto_graf['dist_Gtw']
usar_up = punto_graf['usar_up']
usar_down = punto_graf['usar_down']

# Eje x en log10()
xi = np.array(dist)

# enlace_up
media_up = punto_graf[medida+'_up']
std_up   = punto_graf[medida+'_up'+'_std']
media_up_techo = media_up + std_up
media_up_piso  = media_up - std_up

# ecuacion Enlace up
yi_up = np.array(media_up)
eq_up = girni.linealiza_lstsq(xi[usar_up==1],
                              yi_up[usar_up==1],
                              digitos = precision)
alpha_up = eq_up['alpha']
beta_up  = eq_up['beta']
fd_up = lambda d: -10*alpha_up*(np.log10(d))+beta_up

# Errores up respecto a ecuacion
fi_up = fd_up(xi)
dyi0std  = eq_up['error_std']
punto_graf['fi_up'] = fi_up

# enlace_down
media_down = punto_graf[medida+'_down']
std_down   = punto_graf[medida+'_down'+'_std']
media_down_techo = media_down + std_down
media_down_piso  = media_down - std_down

# ecuación Enlace down
yi_down = np.array(media_down)
eq_down = girni.linealiza_lstsq(xi[usar_down==1],
                                yi_down[usar_down==1],
                                digitos = precision)
alpha_down = eq_down['alpha']
beta_down  = eq_down['beta']
fd_down = lambda d: -10*alpha_down*(np.log10(d))+beta_down

# Errores down respecto a ecuacion
fi_down = fd_down(xi)
dyi1std = eq_down['error_std']
punto_graf['fi_down'] = fi_down

resultado = {'tabla'        : punto_graf,
             'ecuacion_up'  : eq_up,
             'ecuacion_down': eq_down}

ecuacion = resultado.copy()
ecuacion.pop('tabla')
ecuacion = pd.DataFrame(ecuacion)

# SALIDA
pd.options.display.float_format = '{:,.2f}'.format
print(resultado['tabla'])
print('\n ecuaciones: ')
print(ecuacion)
if len(puntoSinDist)>0:
    print('\n Error: Puntos que no registran distancia en coordenadas:')
    print(puntoSinDist)

# graba la tabla resumen rssi
unresumen = carpeta_rsm+'/'+medida+'_'+arch_nombre+'.csv'
resultado['tabla'].to_csv(unresumen)

def graf_gradiente(arch_nombre,punto_graf,ecuacion,var_gen):
    ''' grafica de un gradiente
    '''
    escalabase = 10    # 10
    escala = 'log'
    medida = var_gen['medida']
    precision = var_gen['precision']

    # enlace_up
    media_up = punto_graf[medida+'_up']
    std_up   = punto_graf[medida+'_up'+'_std']
    media_up_techo = media_up + std_up
    media_up_piso  = media_up - std_up
    yi_up = np.array(media_up)
    fi_up = punto_graf['fi_up']

    # enlace_down
    media_down = punto_graf[medida+'_down']
    std_down   = punto_graf[medida+'_down'+'_std']
    media_down_techo = media_down + std_down
    media_down_piso  = media_down - std_down
    yi_down = np.array(media_down)
    fi_down = punto_graf['fi_down']

    dist = punto_graf['dist_Gtw']
    usar_up = punto_graf['usar_up']
    usar_down = punto_graf['usar_down']

    eq_up = ecuacion['ecuacion_up']
    eq_down = ecuacion['ecuacion_down']

    dyi0std = eq_up['error_std']
    dyi1std = eq_down['error_std']

    fig_gradnt,(graf_up,graf_down) = plt.subplots(2,1)
    if escala == 'log':
        graf_up.set_xscale(escala,base=escalabase)
        graf_down.set_xscale(escala,base=escalabase)

    # medida up +/- std
    graf_up.scatter(dist,media_up,color='blue',marker='o')
    graf_up.scatter(dist,media_up_techo,color='blue',marker='+')
    graf_up.scatter(dist,media_up_piso,color='blue',marker='+')

    # medida down +/- std
    graf_down.scatter(dist,media_down,color='orange',marker='o')
    graf_down.scatter(dist,media_down_techo,color='orange',marker='+')
    graf_down.scatter(dist,media_down_piso,color='orange',marker='+')

    # linealizado up
    graf_up.plot(dist,fi_up,label = '  up:'+eq_up['eq_latex'],
                 color='blue')
    graf_up.plot(dist,fi_up+dyi0std,
                 color='blue',linestyle='dotted')
    graf_up.plot(dist,fi_up-dyi0std,
                 color='blue',linestyle='dotted')
    # linealizado down
    graf_down.plot(dist,fi_down,label = 'down:'+eq_down['eq_latex'],
                   color='orange')
    graf_down.plot(dist,fi_down+dyi1std,
                   color='orange',linestyle='dotted')
    graf_down.plot(dist,fi_down-dyi1std,
                   color='orange',linestyle='dotted')

    # ajuste de eje y
    y_min = np.min([np.min(media_up_piso),np.min(media_down_piso)])
    y_max = np.max([np.max(fi_up+dyi0std),np.max(fi_down+dyi1std)])
    if y_min<-135:
        y_min = np.min([np.min([media_up]),np.min([media_down])])
    graf_up.set_ylim(y_min,y_max)
    graf_down.set_ylim(y_min,y_max)

    # etiquetas up
    for unpunto in punto_graf.index:
        media_etiq = np.round(punto_graf[medida+'_up'][unpunto],precision)
        usarpunto  = punto_graf['usar_up'][unpunto]
        
        etiqueta = unpunto #+str(media_etiq)
        xi = punto_graf['dist_Gtw'][unpunto]
        yi = punto_graf['rssi_up'][unpunto]
        if usarpunto:
            graf_up.annotate(etiqueta,(xi,yi), rotation=45)
        if not(usarpunto):
            graf_up.annotate(etiqueta,(xi,yi), rotation=45,color='red')
            graf_up.scatter(dist[unpunto],media_up[unpunto],color='red',marker='o')
            graf_up.scatter(dist[unpunto],media_up_techo[unpunto],
                            color='red',marker='+')

    # etiquetas down
    for unpunto in punto_graf.index:
        media_etiq = np.round(punto_graf[medida+'_down'][unpunto],precision)
        usarpunto = punto_graf['usar_down'][unpunto]
        etiqueta = unpunto #+str(media_etiq)
        xi = punto_graf['dist_Gtw'][unpunto]
        yi = punto_graf['rssi_down'][unpunto]
        if usarpunto:
            graf_down.annotate(etiqueta,(xi,yi), rotation=45)
        if not(usarpunto):
            graf_down.annotate(etiqueta,(xi,yi), rotation=45,color='red')
            graf_down.scatter(dist[unpunto],media_down[unpunto],color='red',marker='o')
            graf_down.scatter(dist[unpunto],media_down_techo[unpunto],
                              color='red',marker='+')

    graf_up.set_ylabel(medida+'_up (dBm)',color='blue')
    graf_up.set_title(arch_nombre+' '+medida)
    graf_up.legend()
    graf_up.grid(True,linestyle='dotted',
                 axis='x', which='both')

    graf_down.set_xlabel('distancia')
    graf_down.set_ylabel(medida+'_down (dBm)',color='brown')
    graf_down.legend()
    graf_down.grid(True,linestyle='dotted',
                   axis='x', which='both')
    return(fig_gradnt)


def graf_diferencia_updown(arch_nombre,punto_graf,var_gen):
    '''Grafica de diferencias up-down
    '''
    medida = var_gen['medida']
    medida_unidad = var_gen['medida_unidad']
    movAvg_cual = var_gen['movAvg_cual']
    movAvg_color = var_gen['movAvg_color']
    precision = var_gen['precision']

    # analiza diferencias de medias y std --------------
    punto_graf['usar_dif'] = punto_graf['usar_up']*punto_graf['usar_down']
    usar_dif = punto_graf['usar_dif']
    punto_graf['dif_mean'] = punto_graf['rssi_up']-punto_graf['rssi_down']
    punto_graf['dif_std']  = punto_graf[medida+'_up_std']-punto_graf[medida+'_down_std']
    punto_graf['dif_top']  = punto_graf['dif_mean']+punto_graf['dif_std']
    punto_graf['dif_bottom'] = punto_graf['dif_mean']-punto_graf['dif_std']

    dist = punto_graf['dist_Gtw']
    usar_dif = punto_graf['usar_dif']

    # linealiza diferencia
    xi_dif = np.array(punto_graf['dist_Gtw'][usar_dif==1])
    yi_dif = np.array(punto_graf['dif_mean'][usar_dif==1])
    eq_dif = girni.linealiza_lstsq(xi_dif,yi_dif,digitos = precision)
    alpha_dif = eq_dif['alpha']
    beta_dif  = eq_dif['beta']
    fdistdif0 = lambda d: -10*alpha_dif*(np.log10(d))+beta_dif
    yi_dif0   = fdistdif0(xi_dif)

    ecua_dif = pd.DataFrame(eq_dif)

    # medias moviles en movAvg_cual[]
    serie_media  = pd.Series(punto_graf['dif_mean'][usar_dif==1])
    movAvg_dmedia = []
    m = len(movAvg_cual)
    for j in range(0,m,1):
        k = movAvg_cual[j]
        movAvg_dmedia.append(serie_media.rolling(k).mean())

    # Grafica diferencia
    escalabase = 10    # 10
    escala = 'log'
    
    fig_dif_updown,graf_dmean = plt.subplots()
    if escala == 'log':
        graf_dmean.set_xscale(escala,base=escalabase)

    # diferencia medida up - down
    graf_dmean.scatter(dist,punto_graf['dif_mean'],color='purple',marker='o')
    graf_dmean.plot(xi_dif,yi_dif0,label = '  dif:'+eq_dif['eq_latex'],color='blue')
    graf_dmean.set_ylabel('media Up-Down (dBm)')

    # medida_up, medias móviles
    for j in range(0,m,1):
        k = str(movAvg_cual[j])
        graf_dmean.plot(dist[usar_dif==1],movAvg_dmedia[j],
                      label='movAvg_dmedia_'+k,color=movAvg_color[j])

    # etiquetas
    i=0
    for unpunto in punto_graf.index:
        media_etiq = np.round(punto_graf[medida+'_up'][unpunto],precision)
        etiqueta = unpunto # + ' ' +str(media_etiq)
        usar_dif = punto_graf['usar_dif'][unpunto]
        xi = punto_graf['dist_Gtw'][unpunto]
        yi = punto_graf.iloc[i]['dif_mean']
        i = i+1
        if usar_dif:
            graf_dmean.annotate(etiqueta,(xi,yi), rotation=30)
        if not(usar_dif):
            graf_dmean.annotate(etiqueta,(xi,yi),color='red',rotation=30)
            graf_dmean.scatter(xi,yi,color='red',marker='o')

    graf_dmean.axhline(0,color='grey')
    graf_dmean.grid(True,linestyle='dotted',axis='x', which='both')

    graf_dmean.legend()
    graf_dmean.set_title(arch_nombre+': '+medida+' media Up-Down')
    return(fig_dif_updown)

def graf_gradiente_std(arch_nombre,punto_graf,medida,precision):
    ''' Grafica de std_up y std_down ------
    '''
    xi_std      = np.array(punto_graf['dist_Gtw'])
    yi_std_up   = np.array(punto_graf[medida+'_up_std'])
    yi_std_down = np.array(punto_graf[medida+'_down_std'])

    # linealiza std
    usar_up   = punto_graf['usar_up']
    usar_down = punto_graf['usar_down']

    eq_std_up = girni.linealiza_lstsq(xi_std[usar_up==1],
                                      yi_std_up[usar_up==1],
                                      digitos = precision)
    alpha_std_up = eq_std_up['alpha']
    beta_std_up  = eq_std_up['beta']
    fdiststd_up0 = lambda d: -10*alpha_std_up*(np.log10(d))+beta_std_up
    yi_std_up0   = fdiststd_up0(xi_std)

    eq_std_down = girni.linealiza_lstsq(xi_std[usar_down==1],
                                        yi_std_down[usar_down==1],
                                        digitos = precision)
    alpha_std_down = eq_std_down['alpha']
    beta_std_down  = eq_std_down['beta']
    fdiststd_down0 = lambda d: -10*alpha_std_down*(np.log10(d))+beta_std_down
    yi_std_down0   = fdiststd_down0(xi_std)

    # grafica std
    escalabase = 10    # 10
    escala = 'log'
    
    fig_grad_std,(graf_std_up,graf_std_down) = plt.subplots(2,1)
    if escala == 'log':
        graf_std_up.set_xscale(escala,base=escalabase)
        graf_std_down.set_xscale(escala,base=escalabase)

    graf_std_up.plot(xi_std,yi_std_up0,
                     label = '  std_up:'+eq_std_up['eq_latex'],
                     color='blue')
        
    # medida up +/- std
    graf_std_up.scatter(xi_std,yi_std_up,color='blue',marker='+')
    graf_std_down.scatter(xi_std,yi_std_down,color='orange',marker='+')

    # etiquetas
    i=0
    for unpunto in punto_graf.index:
        media_etiq = np.round(punto_graf[medida+'_up'][unpunto],precision)
        etiqueta   = unpunto # + ' ' +str(media_etiq)
        usar_up    = punto_graf['usar_up'][unpunto]
        usar_down  = punto_graf['usar_down'][unpunto]
        xi  = punto_graf['dist_Gtw'][unpunto]
        yiu = punto_graf.iloc[i][medida+'_up_std']
        yid = punto_graf.iloc[i][medida+'_down_std']
        i = i+1
        if usar_up:
            graf_std_up.annotate(etiqueta,(xi,yiu), rotation=30)
        if not(usar_up):
            graf_std_up.annotate(etiqueta,(xi,yiu),color='red',rotation=30)
            graf_std_up.scatter(xi,yiu,color='red',marker='o')
        if usar_down:
            graf_std_down.annotate(etiqueta,(xi,yid), rotation=30)
        if not(usar_down):
            graf_std_down.annotate(etiqueta,(xi,yid),color='red',rotation=30)
            graf_std_down.scatter(xi,yid,color='red',marker='o')

    graf_std_down.plot(xi_std,yi_std_down0,
                    label = 'std_down:'+eq_std_down['eq_latex'],
                    color='orange')

    # otras etiquetas
    graf_std_up.set_ylabel('std_up (dBm)', color='blue')
    graf_std_up.set_xlabel('distancia')
    graf_std_up.axhline(0,color='grey')
    graf_std_up.grid(True,linestyle='dotted',
                    axis='x', which='both')
    graf_std_up.set_title(arch_nombre+': '+medida+' std')
    graf_std_up.legend()

    graf_std_down.set_ylabel('std_down (dBm)', color='brown')
    graf_std_down.set_xlabel('distancia')
    graf_std_down.axhline(0,color='grey')
    graf_std_down.grid(True,linestyle='dotted',
                    axis='x', which='both')
    graf_std_down.legend()

    return(fig_grad_std)

fig_gradnt = graf_gradiente(arch_nombre,punto_graf,ecuacion,var_gen)
unarchivo = carpeta_rsm+'/ecuacion_'+arch_nombre+'.png'
fig_gradnt.savefig(unarchivo)

fig_dif_updown = graf_diferencia_updown(arch_nombre,punto_graf,var_gen)
fig_grad_std = graf_gradiente_std(arch_nombre,punto_graf,medida,precision)

plt.show() # pasa a programa principal

4. LoRaWan – Funciones girni_lorawan_lib

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_lorawan_libreria

Instrucciones Python

# Girni LoRaWan librerias 2021-10-28
# LoRaWan, lecturas de Rssi y SNR
# Propuesta: edelros@espol.edu.ec
# http://blog.espol.edu.ec/girni/
import numpy as np
import json
import pandas as pd
import datetime as dt
import os
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import matplotlib.units as munits

import geopandas as gpd
import fiona
import utm

def tablapunto(unarchivoDB,unsensor,gatewayDB,fechainicio,fechafin,zonaGMT=0):
    '''Lectura de archivo.json hacia una tabla0
       selecciona registros de "un sensor"
       graba unreporte.csv con pandas
    '''
    campos = ['domain','entity_id','state','attributes','created']

    # Lectura de archivo json, cambia a formato DataFrame
    tabla0 = pd.DataFrame()
    archivoexiste = os.path.exists(unarchivoDB)
    if archivoexiste:
        tabla0 = pd.read_json(unarchivoDB)
        tabla0 = pd.DataFrame(tabla0,columns=campos)
    else:
        print(' ERROR: NO se encuentra el archivo...')
        print('        revise el nombre de archivo. ')

    # Revisa cada registro 
    leidos = 0
    tabla  = {}
    if archivoexiste:
        # Intervalo de fecha como datetime
        hora_desplaza  = dt.timedelta(hours = abs(zonaGMT))
        fechaformatoDB = '%Y-%m-%dT%H:%M:%S.%f'
        fechaformato   = '%Y-%m-%d %H:%M:%S.%f'
        fechatxt     = fechainicio[0:np.min([26,len(fechainicio)])]
        fechainicio  = dt.datetime.strptime(fechatxt,fechaformato)
        fechatxt     = fechafin[0:np.min([26,len(fechafin)])]
        fechafin     = dt.datetime.strptime(fechatxt,fechaformato)

        # datos de trama
        for cadaregistro in tabla0.index:
            unatrama   = tabla0['attributes'][cadaregistro]
            trama_mqtt = json.loads(unatrama)
            
            # selecciona registro por sensor, en fecha y con datos
            cualsensor = (tabla0['entity_id'][cadaregistro] == unsensor)
            
            unafecha = tabla0['created'][cadaregistro]
            unafecha = unafecha[0:np.min([26,len(unafecha)])]
            unafecha = dt.datetime.strptime(unafecha,fechaformato)
            unafecha = unafecha - hora_desplaza
            
            enfecha  = (unafecha>=fechainicio) and (unafecha<=fechafin)
            condatos = 'applicationID' in trama_mqtt.keys()
            
            selecciona = cualsensor and condatos and enfecha
            if (selecciona == True):        
                # datos en la mensaje MQTT
                publishedAt = trama_mqtt["publishedAt"]
                publishedAt = publishedAt[0:np.min([26,len(publishedAt)])]
                publishedAt = dt.datetime.strptime(publishedAt,fechaformatoDB)
                publishedAt = publishedAt - hora_desplaza
                
                friendly_name = trama_mqtt["friendly_name"]
                deviceName    = trama_mqtt["deviceName"]
                dispositivo   = friendly_name.split('_')[2]
                
                dr    = trama_mqtt["dr"]
                fCnt  = trama_mqtt["fCnt"]
                fPort = trama_mqtt["fPort"]

                objectJSON  = trama_mqtt["objectJSON"]
                datosensor  = json.loads(objectJSON)
                
                freq_tx   = trama_mqtt["txInfo"]["frequency"]
                bandwidth   = trama_mqtt["txInfo"]["loRaModulationInfo"]["bandwidth"]
                spreadingFactor = trama_mqtt["txInfo"]["loRaModulationInfo"]["spreadingFactor"]
                     
                datostrama = {"publishedAt": publishedAt,
                              "dispositivo": dispositivo,
                              "fCnt": fCnt}
                
                # revisa gateway en la red LoRaWan
                tamano = len(trama_mqtt["rxInfo"])
                i = 0
                while i<tamano:
                    gatewayID = trama_mqtt["rxInfo"][i]["gatewayID"]
                    rssi      = trama_mqtt["rxInfo"][i]["rssi"]
                    loRaSNR   = trama_mqtt["rxInfo"][i]["loRaSNR"]
                    channel   = trama_mqtt["rxInfo"][i]["channel"]
                    rfChain   = trama_mqtt["rxInfo"][i]["rfChain"]
                    gtwNum    = gatewayDB[gatewayID]
                    
                    # registra en tabla, incluyendo tramas repetidas
                    datostrama['rssi_up']    = rssi
                    datostrama['snr_up']     = loRaSNR
                    datostrama['channel_up'] = channel
                    datostrama['rfChain_up'] = rfChain
                    datostrama['gtw_rx']     = gtwNum
                    i = i + 1

                # dato del sensor
                equivale = {'Down_datarate': 'dr_down',
                            'Down_rssi'    : 'rssi_down',
                            'Down_snr'     : 'snr_down',
                            'bateria_V'    : 'bateria_V'}
                for undato in datosensor:
                    if undato in equivale.keys():
                        datostrama[equivale[undato]] = datosensor[undato]
                    else:
                        datostrama[undato] = datosensor[undato]

                # datos restantes
                datostrama["frequency"] = freq_tx
                datostrama["bandwidth"] = bandwidth
                datostrama["spreadingFactor"] = spreadingFactor
                datostrama["fPort"]     = fPort
                datostrama["dr"]        = dr
                datostrama["created"]   = unafecha

                leidos = leidos + 1
                
                # revisa registro repetido
                repetido = False
                if leidos>1:
                    trama_num  = (fCnt == tabla[leidos-1]['fCnt'])
                    fecha_pub  = (publishedAt == tabla[leidos-1]['publishedAt'])
                    gtwNum_rep = (gtwNum == tabla[leidos-1]['gtw_rx'])
                    repetido   = (trama_num and fecha_pub and gtwNum_rep)
                if not(repetido):
                    tabla[leidos] = datostrama
                else:
                    leidos = leidos - 1
                
        # convierte diccionario en DataFrame
        if len(tabla)>0:
            tabla = pd.DataFrame(tabla)
            tabla = tabla.T
    return(tabla)

def describe_punto(tabla,codigopunto,carpeta_rsm,variables):
    ''' Descriptores estadisticos de los datos.csv
        de un dispositivo, revisa media, desviació estándar, pmf
        graba unreporte.csv con pandas y genera gráficas
    '''
    medida = variables['medida']
    medida_unidad = variables['medida_unidad']
    medida_normal = variables['medida_normal']
    movAvg_cual = variables['movAvg_cual']
    medida_grafica = variables['medida_grafica']
    movAvg_color = variables['movAvg_color']
    guarda = variables['guarda']
    precision = variables['precision']
                   
    fechaformato = "%Y-%m-%d %H:%M:%S.%f"
    # medida intervalo
    medida_min = np.min(medida_normal)
    medida_max = np.max(medida_normal)

    # fechas series a datetime
    tabla['publishedAt'] = pd.to_datetime(tabla['publishedAt'],
                                          format=fechaformato)
    tabla['created'] = pd.to_datetime(tabla['publishedAt'],
                                      format=fechaformato)

    # revisa errores de medida
    tabla["error_up"]   = 0
    tabla["error_down"] = 0
    for undato in tabla['publishedAt'].keys():   
        medida_up = tabla[medida+'_up'][undato]
        enrango = (medida_up>=np.min(medida_normal))
        enrango = (enrango and medida_up<=np.max(medida_normal))
        if not(enrango):
            tabla.at[undato,"error_up"] = 1
        
        medida_down = tabla[medida+'_down'][undato]
        enrango = (medida_down>=np.min(medida_normal))
        enrango = (enrango and medida_down<=np.max(medida_normal))
        if not(enrango):
            tabla.at[undato,"error_down"] = 1      

    # tasa error trama
    leidos = len(tabla)
    if leidos > 0:
        error_up   = np.sum(tabla['error_up'])
        error_up   = error_up/leidos
        error_down = np.sum(tabla['error_down'])
        error_down = error_down/leidos

    # descriptor estadístico, datos sin errores
    condicion_up = (tabla['error_up']==0)
    condicion_down = (tabla['error_down']==0)

    medida_up = tabla[condicion_up][medida+'_up']
    describe_up = medida_up.describe()
    describe_up['error_trama'] = error_up

    medida_down = tabla[condicion_down][medida+'_down']
    describe_down = medida_down.describe()
    describe_down['error_trama'] = error_down

    descriptor = describe_up.copy()
    descriptor = pd.concat([descriptor,describe_down],axis=1)
    descriptor['dispositivo'] = tabla['dispositivo'][0]

    pmf_up   = medida_pmf(medida_up,describe_up)
    pmf_down = medida_pmf(medida_down,describe_down)

    pmf_punto = {'pmf':{'pmf_up'   : pmf_up,
                        'pmf_down' : pmf_down}}
    pmf_punto = pd.DataFrame(pmf_punto)
    pmf_punto['dispositivo'] = tabla['dispositivo'][0]

    # Para gráficas
    # medias moviles en movAvg_cual[]
    serie_up  = pd.Series(medida_up)
    movAvg_up_mean = []
    movAvg_up_std = []
    m = len(movAvg_cual)
    for j in range(0,m,1):
        k = movAvg_cual[j]
        movAvg_up_mean.append(list(serie_up.rolling(k).mean()))
        movAvg_up_std.append(list(serie_up.rolling(k).std()))
        
    serie_down = pd.Series(medida_down)
    movAvg_down_mean = []
    movAvg_down_std = []
    for j in range(0,m,1):
        k = movAvg_cual[j]
        movAvg_down_mean.append(list(serie_down.rolling(k).mean()))
        movAvg_down_std.append(list(serie_down.rolling(k).std()))

    movAvgData ={'movAvg_cual'   : movAvg_cual,
                 'movAvg_color'  : movAvg_color,
                 'movAvg_up_mean'  : movAvg_up_mean,
                 'movAvg_down_mean': movAvg_down_mean,
                 'movAvg_up_std'   : movAvg_up_std,
                 'movAvg_down_std' : movAvg_down_std
                 }

    grafData ={'codigopunto' : codigopunto,
               'medida' : medida,
               'precision': precision,
               'medida_unidad' : medida_unidad,
               'medida_grafica': medida_grafica
               }
    resultado = [tabla,descriptor,pmf_punto,movAvgData,grafData]
    return(resultado)

# función de probabilidad de masa pmf
def medida_pmf(valores,undescriptor):
    pmin   = np.min(valores)
    pmax   = np.max(valores)
    tramo  = int(pmax-pmin)
    conteo = np.zeros(tramo+1,dtype=int)
    intervalo = np.arange(pmin,pmax+1,1)
    for valor in valores:
        donde = np.where(intervalo == valor)
        conteo[donde] = conteo[donde] + 1
    freq_relativa = np.array(conteo)/np.sum(conteo)
    unpmf = {'intervalo' : list(intervalo),
             'freq_relativa' : list(freq_relativa)}
    return(unpmf)

# GRAFICA -----
def graf_puntos_serie(tabla,descriptor,movAvgData,grafData):
    ''' grafica la serie de tiempo de cada punto
        añade medias móviles
    '''
    # ajuste de formato de fecha para eje x
    converter = mdates.ConciseDateConverter()
    munits.registry[np.datetime64] = converter
    munits.registry[dt.date] = converter
    munits.registry[dt.datetime] = converter

    # datos para grafica
    precision   = grafData['precision']
    medida = grafData['medida']
    medida_unidad  = grafData['medida_unidad']
    medida_grafica = grafData['medida_grafica']
    
    movAvg_cual  = movAvgData['movAvg_cual']
    movAvg_color = movAvgData['movAvg_color']
    
    media_up    = descriptor[medida+'_up']['mean']
    std_up      = descriptor[medida+'_up']['std']
    media_down  = descriptor[medida+'_down']['mean']
    std_down    = descriptor[medida+'_down']['std']

    # ajuste de intervalo eje y
    y_min = np.min([np.min(medida_grafica),
                    media_up - 2*std_up,
                    media_down - 2*std_down])
    y_max = np.max([np.max(medida_grafica),
                    media_up + 2*std_up,
                    media_down + 2*std_down])
    
    # selecciona sin error
    condicion_up = (tabla['error_up']==0)
    condicion_down = (tabla['error_down']==0)

    # grafica
    fig_serie,(graf_up,graf_down) = plt.subplots(2,1)
    
    # medida_up -----
    graf_up.plot(tabla[condicion_up]['publishedAt'],
                 tabla[condicion_up][medida+'_up'],
                 color='blue',marker ='.',
                 linestyle='')
    
    # medida_up, medias y std
    etiq_up = str(np.round(media_up,precision))+' +/- '
    etiq_up = etiq_up + str(np.round(std_up,precision))
    graf_up.axhline(media_up,
                    color='blue',label=etiq_up)
    graf_up.axhline(media_up-std_up,
                    color='blue',linestyle='dotted')
    graf_up.axhline(media_up+std_up,
                    color='blue',linestyle='dotted')
    
    # medida_up, medias móviles
    m = len(movAvg_cual)
    for j in range(0,m,1):
        k = str(movAvg_cual[j])
        graf_up.plot(tabla[condicion_up]['publishedAt'],
                     movAvgData['movAvg_up_mean'][j],
                     label='movAvg_'+k,
                     color=movAvg_color[j])
    
    graf_up.set_ylim(y_min,y_max)
    graf_up.set_ylabel(medida+'_up ('+medida_unidad+')',
                       color='blue')
    graf_up.legend()
    graf_up.grid(True,linestyle='dotted',
                 axis='x',which='both')

    # medida_down -------
    graf_down.plot(tabla[condicion_down]['publishedAt'],
                   tabla[condicion_down][medida+'_down'],
                   color='brown',marker ='.',
                   linestyle='')

    # medida_down, medias y std
    etiq_down = str(np.round(media_down,precision))+' +/- '
    etiq_down = etiq_down + str(np.round(std_down,precision))
    graf_down.axhline(media_down,
                      color='brown',label=etiq_down)
    graf_down.axhline(media_down+std_down,
                      color='brown',linestyle='dotted')
    graf_down.axhline(media_down-std_down,
                      color='brown',linestyle='dotted')
    
    # medida_down, medias moviles
    for j in range(0,m,1):
        k = str(movAvg_cual[j])
        graf_down.plot(tabla[condicion_down]['publishedAt'],
                       movAvgData['movAvg_down_mean'][j],
                       label='movAvg_'+k,
                       color=movAvg_color[j])
    
    graf_down.set_ylim(y_min,y_max)
    graf_down.set_xlabel('fecha')
    graf_down.set_ylabel(medida+'_down ('+medida_unidad+')',
                         color='brown')
    graf_down.legend()
    graf_down.grid(True,linestyle='dotted',
                   axis='x', which='both')
    graf_up.set_title('Serie: '+grafData['codigopunto']+' '+ medida)
    plt.tight_layout()
    return(fig_serie)

def graf_puntos_pmf(pmf_punto,descriptor,grafData):
    ''' grafica función de probabilida de masa
        para cada punto, media +/- std
    '''
    # datos para grafica
    x_pmfup   = pmf_punto['pmf']['pmf_up']['intervalo']
    y_pmfup   = pmf_punto['pmf']['pmf_up']['freq_relativa']
    x_pmfdown = pmf_punto['pmf']['pmf_down']['intervalo']
    y_pmfdown = pmf_punto['pmf']['pmf_down']['freq_relativa']
    
    precision = grafData['precision']
    medida = grafData['medida']
    medida_unidad  = grafData['medida_unidad']
    medida_grafica = grafData['medida_grafica']
    
    media_up   = descriptor[medida+'_up']['mean']
    std_up     = descriptor[medida+'_up']['std']
    media_down = descriptor[medida+'_down']['mean']
    std_down   = descriptor[medida+'_down']['std']

    prob_max = 0.40
    # ajuste de intervalo eje y
    y_min = np.min([np.min(medida_grafica),
                    media_up - 2*std_up,
                    media_down - 2*std_down])
    y_max = np.max([np.max(medida_grafica),
                    media_up + 2*std_up,
                    media_down + 2*std_down])
    # grafica
    fig_pmf,graf_pmf = plt.subplots()
    etiq_up = str(np.round(media_up,precision)) +' +/- '
    etiq_up = etiq_up + str(np.round(std_up,precision))
    graf_pmf.plot(x_pmfup,y_pmfup,
                  label='media_up '+etiq_up,
                  color='blue')
    graf_pmf.axvline(media_up,color='blue')
    graf_pmf.axvline(media_up+std_up,
                     linestyle='dotted',color='blue')
    graf_pmf.axvline(media_up-std_up,
                     linestyle='dotted',color='blue')

    etiq_down = str(np.round(media_down,precision))+' +/- '
    etiq_down = etiq_down + str(np.round(std_down,precision))
    graf_pmf.plot(x_pmfdown,y_pmfdown,
                  label='media_down '+etiq_down,
                  color='brown')
    graf_pmf.axvline(media_down,color='brown')
    graf_pmf.axvline(media_down+std_down,
                     linestyle='dotted',color='brown')
    graf_pmf.axvline(media_down-std_down,
                     linestyle='dotted',color='brown')

    graf_pmf.set_title('pmf: '+grafData['codigopunto']+' '+medida)
    graf_pmf.set_xlim(y_min,y_max)
    graf_pmf.set_ylim(0,prob_max)
    graf_pmf.set_xlabel(medida+' ('+medida_unidad+')')
    graf_pmf.set_ylabel('frecuencia relativa')
    graf_pmf.legend()
    graf_pmf.grid(True,linestyle='dotted',
                  axis='x', which='both')
    return(fig_pmf)

def graf_puntos_std(tabla,descriptor,movAvgData,grafData):
    ''' grafica serie de std usando medias moviles
        para cada punto, media_std
    '''
    # ajuste de formato de fecha para eje x
    converter = mdates.ConciseDateConverter()
    munits.registry[np.datetime64] = converter
    munits.registry[dt.date] = converter
    munits.registry[dt.datetime] = converter
    
    # datos para grafica
    precision = grafData['precision']
    medida    = grafData['medida']
    medida_unidad = grafData['medida_unidad']

    movAvg_cual = movAvgData['movAvg_cual']
    movAvg_color = movAvgData['movAvg_color']

    # selecciona sin error
    condicion_up   = (tabla['error_up']==0)
    condicion_down = (tabla['error_down']==0)
    
    # ajuste de intervalo eje y
    y_min = 0
    y_max = np.max([2, 2*descriptor[medida+'_up']['std'],
                    2*descriptor[medida+'_down']['std']])
    # grafica
    fig_std,(graf_stdUp,graf_stdDown) = plt.subplots(2,1)
    
    # std up
    std_up = np.round(descriptor[medida+'_up']['std'],precision)
    graf_stdUp.axhline(std_up,label='std '+str(std_up),
                       color='blue')
    m = len(movAvg_cual)
    for j in range(0,m,1):
        k = str(movAvg_cual[j])
        graf_stdUp.plot(tabla[condicion_up]['publishedAt'],
                        movAvgData['movAvg_up_std'][j],
                        label='movAvg_'+k,
                        color=movAvg_color[j])
    graf_stdUp.set_ylim(y_min,y_max)
    graf_stdUp.set_ylabel('std_up ('+medida_unidad+')',
                          color='blue')
    graf_stdUp.legend()
    graf_stdUp.grid(True,linestyle='dotted',
                    axis='x', which='both')
    graf_stdUp.set_title('std: '+grafData['codigopunto']+' '+ medida)

    # std down
    std_down = np.round(descriptor[medida+'_down']['std'],precision)
    graf_stdDown.axhline(std_down,label='std '+str(std_down),
                         color='brown')
    for j in range(0,m,1):
        k = str(movAvg_cual[j])
        graf_stdDown.plot(tabla[condicion_down]['publishedAt'],
                          movAvgData['movAvg_down_std'][j],
                          label='movAvg_'+k,color=movAvg_color[j])
    graf_stdDown.set_ylim(y_min,y_max)
    graf_stdDown.set_xlabel('fecha')
    graf_stdDown.set_ylabel('std_down ('+medida_unidad+')',
                            color='brown')
    graf_stdDown.legend()
    graf_stdDown.grid(True,linestyle='dotted',
                      axis='x', which='both')
    plt.tight_layout()
    return(fig_std)


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

    # coeficiente de correlación
    coeficientes = np.corrcoef(xilog,yi)
    r = coeficientes[0,1]
    # coeficiente de determinación
    r2 = r**2

    # 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)
    
    yimin = np.around(np.min(yi),2)
    yimax = np.around(np.max(yi),2)
    
    unaecuacion = {'alpha'   : alpha,
                   'beta'    : beta,
                   'coef_correlacion'   : r,
                   'coef_determinacion' : r2,
                   'eq_latex': fdtxt0,
                   'intervalox' : [np.min(xi),np.max(xi)],
                   'error_medio': dyi0mean,
                   'error_std'  : dyi0std,
                   'eqg_latex'  : grtxt0,
                   'intervaloy' : [yimin,yimax],
                   'errorx_medio': dxi0mean,
                   'errorx_std'  : dxi0std,
                   }
    return(unaecuacion)

def ecuacion_gradiente(carpeta_rsm,arch_coord,var_gen):
    ''' desarrolla la ecuación para un gradiente
    '''
    var_gen['movAvg_cual'] = [2,4] #cada cuantas muestras
    partes = carpeta_rsm.strip('/').split('_')
    arch_nombre = partes[1]+'_'+partes[2]

    carp_coord = var_gen['carp_coord']
    medida = var_gen['medida']
    precision = var_gen['precision']
    
    # lee coordenadas y su distancia al gateway
    dist_Gtw = coordenadas_leer(arch_coord,carp_coord)

    # crea lista.txt de archivos a usar si no existe
    arch_lista = medida+"_"+arch_nombre+"_lista.txt"
    archivo_ruta  = carpeta_rsm + '/' + arch_lista
    archivoexiste = os.path.exists(archivo_ruta)
    if not(archivoexiste):
        puntoUsar = pd.DataFrame(columns=['punto', 'up', 'down'])
        for unarchivo in os.listdir(carpeta_rsm):
            verifica = unarchivo.startswith('describe_')
            verifica = verifica and unarchivo.endswith('.csv')
            if verifica:
                partes  = unarchivo.strip('.csv').split('_')
                unpunto = partes[2]
                puntoUsar = puntoUsar.append({'punto' : unpunto,
                                              'up' : 1,
                                              'down' : 1},
                                             ignore_index = True)
        puntoUsar = puntoUsar.set_index('punto')
        puntoUsar.to_csv(archivo_ruta)

    # usa lista.txt de archivos seleccionados 1 ó 0 para enlace up,down
    if (archivoexiste):
        puntoUsar = pd.read_csv(archivo_ruta)
        puntoUsar = puntoUsar.set_index('punto')
      
    # Datos para gráfica desde lista.txt
    punto_graf  = pd.DataFrame()
    puntoSinDist = []
    for unarchivo in os.listdir(carpeta_rsm):
        if unarchivo.startswith('describe_'):
            codigopunto = unarchivo.strip('.csv').split('_')[2]
            
            # lectura del archivo
            unresumen  = carpeta_rsm+'/'+unarchivo
            descriptor = pd.read_csv(unresumen,index_col='Unnamed: 0')

            if (codigopunto in dist_Gtw.index):
                undato = {'codigopunto': codigopunto,
                          'dist_Gtw'   : dist_Gtw[codigopunto],
                          medida+'_up' : descriptor['rssi_up']['mean'],
                          medida+'_up'+'_std': descriptor['rssi_up']['std'],
                          'usar_up' : puntoUsar['up'][codigopunto],
                          medida+'_down': descriptor['rssi_down']['mean'],
                          medida+'_down'+'_std' :descriptor['rssi_down']['std'],
                          'usar_down' : puntoUsar['down'][codigopunto],
                          'dispositivo': descriptor['dispositivo'][0]
                          }
                punto_graf = punto_graf.append(undato,ignore_index=True)
            else:
                puntoSinDist.append(codigopunto)
            
    punto_graf = punto_graf.set_index('codigopunto')
    punto_graf = punto_graf.sort_values('dist_Gtw' )

    # selecciona datos
    dist = punto_graf['dist_Gtw']
    usar_up = punto_graf['usar_up']
    usar_down = punto_graf['usar_down']

    # Eje x en log10()
    xi = np.array(dist)

    # enlace_up
    media_up = punto_graf[medida+'_up']
    std_up   = punto_graf[medida+'_up'+'_std']
    media_up_techo = media_up + std_up
    media_up_piso  = media_up - std_up

    # ecuacion enlace up
    yi_up = np.array(media_up)
    eq_up = linealiza_lstsq(xi[usar_up==1],
                            yi_up[usar_up==1],
                            digitos = precision)
    alpha_up = eq_up['alpha']
    beta_up  = eq_up['beta']
    fd_up = lambda d: -10*alpha_up*(np.log10(d))+beta_up

    # errores up de puntos vs ecuacion
    fi_up = fd_up(xi)
    dyi0std  = eq_up['error_std']
    punto_graf['fi_up'] = fi_up

    # enlace_down
    media_down = punto_graf[medida+'_down']
    std_down   = punto_graf[medida+'_down'+'_std']
    media_down_techo = media_down + std_down
    media_down_piso  = media_down - std_down

    # ecuación Enlace down
    yi_down = np.array(media_down)
    eq_down = linealiza_lstsq(xi[usar_down==1],
                              yi_down[usar_down==1],
                              digitos = precision)
    alpha_down = eq_down['alpha']
    beta_down  = eq_down['beta']
    fd_down = lambda d: -10*alpha_down*(np.log10(d))+beta_down

    # errores down de puntos vs ecuacion
    fi_down = fd_down(xi)
    dyi1std = eq_down['error_std']
    punto_graf['fi_down'] = fi_down

    resultado = {'tabla'         : punto_graf,
                 'ecuacion_up'   : eq_up,
                 'ecuacion_down' : eq_down,
                 'puntoSinDist'  : puntoSinDist}
    
    return (resultado)

def coordenadas_leer(arch_coord,carp_coord):
    ''' lista las coordenadas desde un archivo
        ubicado en la carpeta
    '''
    if len(carp_coord)>0:
        carp_coord = carp_coord+'/'
    archivo_ruta = carp_coord+arch_coord
    arch_coord_existe = os.path.exists(archivo_ruta)

    if arch_coord_existe:
        puntos = pd.read_csv(archivo_ruta,index_col='punto')
        # busca gateway y ubicación de punto observado
        Gw_nombre = ''
        for unpunto in puntos.index:
            if unpunto.startswith('Gw'):
                Gw_nombre = unpunto
        dist_Gtw = puntos['dist_'+Gw_nombre]   
    else:
        print(' ERROR: NO se encuentra el archivo...')
        print('        revise el ruta/nombre de archivo. ')
    return(dist_Gtw)

def graf_gradiente(arch_nombre,punto_graf,ecuacion,var_gen):
    ''' grafica de un gradiente
    '''
    escalabase = 10    # 10
    escala = 'log'
    medida = var_gen['medida']
    precision = var_gen['precision']

    # enlace_up
    media_up = punto_graf[medida+'_up']
    std_up   = punto_graf[medida+'_up'+'_std']
    media_up_techo = media_up + std_up
    media_up_piso  = media_up - std_up
    yi_up = np.array(media_up)
    fi_up = punto_graf['fi_up']

    # enlace_down
    media_down = punto_graf[medida+'_down']
    std_down   = punto_graf[medida+'_down'+'_std']
    media_down_techo = media_down + std_down
    media_down_piso  = media_down - std_down
    yi_down = np.array(media_down)
    fi_down = punto_graf['fi_down']

    dist = punto_graf['dist_Gtw']
    usar_up = punto_graf['usar_up']
    usar_down = punto_graf['usar_down']

    eq_up = ecuacion['ecuacion_up']
    eq_down = ecuacion['ecuacion_down']

    dyi0std = eq_up['error_std']
    dyi1std = eq_down['error_std']

    fig_gradnt,(graf_up,graf_down) = plt.subplots(2,1)
    if escala == 'log':
        graf_up.set_xscale(escala,base=escalabase)
        graf_down.set_xscale(escala,base=escalabase)

    # medida up +/- std
    graf_up.scatter(dist,media_up,color='blue',marker='o')
    graf_up.scatter(dist,media_up_techo,color='blue',marker='+')
    graf_up.scatter(dist,media_up_piso,color='blue',marker='+')

    # medida down +/- std
    graf_down.scatter(dist,media_down,color='orange',marker='o')
    graf_down.scatter(dist,media_down_techo,color='orange',marker='+')
    graf_down.scatter(dist,media_down_piso,color='orange',marker='+')

    # linealizado up
    graf_up.plot(dist,fi_up,label = '  up:'+eq_up['eq_latex'],
                 color='blue')
    graf_up.plot(dist,fi_up+dyi0std,
                 color='blue',linestyle='dotted')
    graf_up.plot(dist,fi_up-dyi0std,
                 color='blue',linestyle='dotted')
    # linealizado down
    graf_down.plot(dist,fi_down,label = 'down:'+eq_down['eq_latex'],
                   color='orange')
    graf_down.plot(dist,fi_down+dyi1std,
                   color='orange',linestyle='dotted')
    graf_down.plot(dist,fi_down-dyi1std,
                   color='orange',linestyle='dotted')

    # ajuste de eje y
    y_min = np.min([np.min(media_up_piso),np.min(media_down_piso)])
    y_max = np.max([np.max(yi_up+dyi0std),np.max(yi_down+dyi1std)])
    if y_min<-135:
        y_min = np.min([np.min([media_up]),np.min([media_down])])
    graf_up.set_ylim(y_min,y_max)
    graf_down.set_ylim(y_min,y_max)

    # etiquetas up
    for unpunto in punto_graf.index:
        media_etiq = np.round(punto_graf[medida+'_up'][unpunto],precision)
        usarpunto  = punto_graf['usar_up'][unpunto]
        
        etiqueta = unpunto #+str(media_etiq)
        xi = punto_graf['dist_Gtw'][unpunto]
        yi = punto_graf['rssi_up'][unpunto]
        if usarpunto:
            graf_up.annotate(etiqueta,(xi,yi), rotation=45)
        if not(usarpunto):
            graf_up.annotate(etiqueta,(xi,yi), rotation=45,color='red')
            graf_up.scatter(dist[unpunto],media_up[unpunto],color='red',marker='o')
            graf_up.scatter(dist[unpunto],media_up_techo[unpunto],
                            color='red',marker='+')

    # etiquetas down
    for unpunto in punto_graf.index:
        media_etiq = np.round(punto_graf[medida+'_down'][unpunto],precision)
        usarpunto = punto_graf['usar_down'][unpunto]
        etiqueta = unpunto #+str(media_etiq)
        xi = punto_graf['dist_Gtw'][unpunto]
        yi = punto_graf['rssi_down'][unpunto]
        if usarpunto:
            graf_down.annotate(etiqueta,(xi,yi), rotation=45)
        if not(usarpunto):
            graf_down.annotate(etiqueta,(xi,yi), rotation=45,color='red')
            graf_down.scatter(dist[unpunto],media_down[unpunto],color='red',marker='o')
            graf_down.scatter(dist[unpunto],media_down_techo[unpunto],
                              color='red',marker='+')

    graf_up.set_ylabel(medida+'_up (dBm)',color='blue')
    graf_up.set_title(arch_nombre+' '+medida)
    graf_up.legend()
    graf_up.grid(True,linestyle='dotted',
                 axis='x', which='both')

    graf_down.set_xlabel('distancia')
    graf_down.set_ylabel(medida+'_down (dBm)',color='brown')
    graf_down.legend()
    graf_down.grid(True,linestyle='dotted',
                   axis='x', which='both')
    return(fig_gradnt)

def graf_diferencia_updown(arch_nombre,punto_graf,var_gen):
    '''Grafica de diferencias up-down
    '''
    medida = var_gen['medida']
    medida_unidad = var_gen['medida_unidad']
    movAvg_cual = var_gen['movAvg_cual']
    movAvg_color = var_gen['movAvg_color']
    precision = var_gen['precision']

    # analiza diferencias de medias y std --------------
    punto_graf['usar_dif'] = punto_graf['usar_up']*punto_graf['usar_down']
    usar_dif = punto_graf['usar_dif']
    punto_graf['dif_mean'] = punto_graf['rssi_up']-punto_graf['rssi_down']
    punto_graf['dif_std']  = punto_graf[medida+'_up_std']-punto_graf[medida+'_down_std']
    punto_graf['dif_top']  = punto_graf['dif_mean']+punto_graf['dif_std']
    punto_graf['dif_bottom'] = punto_graf['dif_mean']-punto_graf['dif_std']

    dist = punto_graf['dist_Gtw']
    usar_dif = punto_graf['usar_dif']

    # linealiza diferencia
    xi_dif = np.array(punto_graf['dist_Gtw'][usar_dif==1])
    yi_dif = np.array(punto_graf['dif_mean'][usar_dif==1])
    eq_dif = linealiza_lstsq(xi_dif,yi_dif,digitos = precision)
    alpha_dif = eq_dif['alpha']
    beta_dif  = eq_dif['beta']
    fdistdif0 = lambda d: -10*alpha_dif*(np.log10(d))+beta_dif
    yi_dif0   = fdistdif0(xi_dif)

    ecua_dif = pd.DataFrame(eq_dif)

    # medias moviles en movAvg_cual[]
    serie_media  = pd.Series(punto_graf['dif_mean'][usar_dif==1])
    movAvg_dmedia = []
    m = len(movAvg_cual)
    for j in range(0,m,1):
        k = movAvg_cual[j]
        movAvg_dmedia.append(serie_media.rolling(k).mean())

    # Grafica diferencia
    escalabase = 10    # 10
    escala = 'log'
    
    fig_dif_updown,graf_dmean = plt.subplots()
    if escala == 'log':
        graf_dmean.set_xscale(escala,base=escalabase)

    # diferencia medida up - down
    graf_dmean.scatter(dist,punto_graf['dif_mean'],color='purple',marker='o')
    graf_dmean.plot(xi_dif,yi_dif0,label = '  dif:'+eq_dif['eq_latex'],color='blue')
    graf_dmean.set_ylabel('media Up-Down (dBm)')

    # medida_up, medias móviles
    for j in range(0,m,1):
        k = str(movAvg_cual[j])
        graf_dmean.plot(dist[usar_dif==1],movAvg_dmedia[j],
                      label='movAvg_dmedia_'+k,color=movAvg_color[j])

    # etiquetas
    i=0
    for unpunto in punto_graf.index:
        media_etiq = np.round(punto_graf[medida+'_up'][unpunto],precision)
        etiqueta = unpunto # + ' ' +str(media_etiq)
        usar_dif = punto_graf['usar_dif'][unpunto]
        xi = punto_graf['dist_Gtw'][unpunto]
        yi = punto_graf.iloc[i]['dif_mean']
        i = i+1
        if usar_dif:
            graf_dmean.annotate(etiqueta,(xi,yi), rotation=30)
        if not(usar_dif):
            graf_dmean.annotate(etiqueta,(xi,yi),color='red',rotation=30)
            graf_dmean.scatter(xi,yi,color='red',marker='o')

    graf_dmean.axhline(0,color='grey')
    graf_dmean.grid(True,linestyle='dotted',axis='x', which='both')

    graf_dmean.legend()
    graf_dmean.set_title(arch_nombre+': '+medida+' media Up-Down')
    return(fig_dif_updown)

def graf_gradiente_std(arch_nombre,punto_graf,medida,precision):
    ''' Grafica de std_up y std_down ------
    '''
    xi_std      = np.array(punto_graf['dist_Gtw'])
    yi_std_up   = np.array(punto_graf[medida+'_up_std'])
    yi_std_down = np.array(punto_graf[medida+'_down_std'])

    # linealiza std
    usar_up   = punto_graf['usar_up']
    usar_down = punto_graf['usar_down']

    eq_std_up = linealiza_lstsq(xi_std[usar_up==1],
                                      yi_std_up[usar_up==1],
                                      digitos = precision)
    alpha_std_up = eq_std_up['alpha']
    beta_std_up  = eq_std_up['beta']
    fdiststd_up0 = lambda d: -10*alpha_std_up*(np.log10(d))+beta_std_up
    yi_std_up0   = fdiststd_up0(xi_std)

    eq_std_down = linealiza_lstsq(xi_std[usar_down==1],
                                        yi_std_down[usar_down==1],
                                        digitos = precision)
    alpha_std_down = eq_std_down['alpha']
    beta_std_down  = eq_std_down['beta']
    fdiststd_down0 = lambda d: -10*alpha_std_down*(np.log10(d))+beta_std_down
    yi_std_down0   = fdiststd_down0(xi_std)

    # grafica std
    escalabase = 10    # 10
    escala = 'log'
    
    fig_grad_std,(graf_std_up,graf_std_down) = plt.subplots(2,1)
    if escala == 'log':
        graf_std_up.set_xscale(escala,base=escalabase)
        graf_std_down.set_xscale(escala,base=escalabase)

    graf_std_up.plot(xi_std,yi_std_up0,
                     label = '  std_up:'+eq_std_up['eq_latex'],
                     color='blue')
        
    # medida up +/- std
    graf_std_up.scatter(xi_std,yi_std_up,color='blue',marker='+')
    graf_std_down.scatter(xi_std,yi_std_down,color='orange',marker='+')

    # etiquetas
    i=0
    for unpunto in punto_graf.index:
        media_etiq = np.round(punto_graf[medida+'_up'][unpunto],precision)
        etiqueta   = unpunto # + ' ' +str(media_etiq)
        usar_up    = punto_graf['usar_up'][unpunto]
        usar_down  = punto_graf['usar_down'][unpunto]
        xi  = punto_graf['dist_Gtw'][unpunto]
        yiu = punto_graf.iloc[i][medida+'_up_std']
        yid = punto_graf.iloc[i][medida+'_down_std']
        i = i+1
        if usar_up:
            graf_std_up.annotate(etiqueta,(xi,yiu), rotation=30)
        if not(usar_up):
            graf_std_up.annotate(etiqueta,(xi,yiu),color='red',rotation=30)
            graf_std_up.scatter(xi,yiu,color='red',marker='o')
        if usar_down:
            graf_std_down.annotate(etiqueta,(xi,yid), rotation=30)
        if not(usar_down):
            graf_std_down.annotate(etiqueta,(xi,yid),color='red',rotation=30)
            graf_std_down.scatter(xi,yid,color='red',marker='o')

    graf_std_down.plot(xi_std,yi_std_down0,
                    label = 'std_down:'+eq_std_down['eq_latex'],
                    color='orange')
    # otras etiquetas
    graf_std_up.set_ylabel('std_up (dBm)', color='blue')
    graf_std_up.set_xlabel('distancia')
    graf_std_up.axhline(0,color='grey')
    graf_std_up.grid(True,linestyle='dotted',
                    axis='x', which='both')
    graf_std_up.set_title(arch_nombre+': '+medida+' std')
    graf_std_up.legend()

    graf_std_down.set_ylabel('std_down (dBm)', color='brown')
    graf_std_down.set_xlabel('distancia')
    graf_std_down.axhline(0,color='grey')
    graf_std_down.grid(True,linestyle='dotted',
                    axis='x', which='both')
    graf_std_down.legend()

    return(fig_grad_std)

def coordenadas_kml_utm(archivo_ruta,zonaNum,zonaLetra,precision):
    # Lectura de archivo
    gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'
    puntos  = gpd.read_file(archivo_ruta, driver='KML')
    puntos['dist'] = 0

    tabla = pd.DataFrame()
    for i in puntos.index:
        # en latitud, longitud y altura
        punto_nombre = puntos.loc[i]['Name']
        longitud = puntos.loc[i]['geometry'].x
        latitud  = puntos.loc[i]['geometry'].y
        altitud  = puntos.loc[i]['geometry'].z
        altitud  = np.round(altitud,precision)
        #en UTM
        coord_utm  = utm.from_latlon(latitud,longitud,
                                     zonaNum,zonaLetra)
        utm_este  = np.round(coord_utm[0],precision)
        utm_norte = np.round(coord_utm[1],precision)
        
        # distancias a Gateways
        dist = 0.0
        if punto_nombre.startswith('Gw'):
            n_Gtw = punto_nombre
            x1 = utm_este
            y1 = utm_norte
        else:
            x2 = utm_este
            y2 = utm_norte
            dist  = np.sqrt((x2-x1)**2 + (y2-y1)**2)
        dist = np.round(dist,precision)
        
        unpunto = {'punto'       : punto_nombre,
                   'utm_este'    : utm_este,
                   'utm_norte'   : utm_norte,
                   'utm_zNum'    : int(zonaNum),
                   'utm_zLetra'  : zonaLetra,
                   'altitud'     : altitud,
                   'altitud_gps' : altitud,
                   'dist_'+n_Gtw : dist,
                   'latitud'     : latitud,
                   'longitud'    : longitud
                   }
        tabla = tabla.append(unpunto,ignore_index=True)
    tabla = tabla.set_index('punto')
    return(tabla)

def revisa_perfil(archivo_ruta,alturaGw,alturapunto,
                  plantas,muestras=41,guarda=False):
    LineaObs = archivo_ruta.strip('.csv').split('_')[1]

    planta_desde  = plantas[0]
    planta_hasta  = plantas[1]
    planta_altura = plantas[2]

    # lectura de datos
    puntos = pd.read_csv(archivo_ruta,index_col='punto')

    # busca gateway
    n_Gtw = ''
    for unpunto in puntos.index:
        if unpunto.startswith('Gw'):
            n_Gtw = unpunto

    # perfil  
    xi = np.array(puntos['dist_'+n_Gtw])
    yi = np.array(puntos['altitud'])
    etiqueta = puntos.index

    # plantaciones
    conplantas = (xi > planta_desde) & (xi < planta_hasta)
    plantacion = yi[conplantas]+planta_altura

    # GRAFICA
    # plantacion
    plt.fill_between(xi[conplantas],
                     yi[conplantas],
                     plantacion,
                     color='lightgreen')
    # perfil
    xi = puntos['dist_'+n_Gtw]
    yi = puntos['altitud']
    plt.plot(xi,yi,color='brown')
        
    # antenas
    for unpunto in puntos.index:
        xip = xi[unpunto]
        yip = yi[unpunto]
        if unpunto.startswith('Gw'):
            xigtw = xip
            yigtw = yip+alturaGw
            plt.scatter(xip ,yigtw)
            plt.annotate(unpunto,(xigtw,yigtw))
        else:
            xid = xip
            yid = yip+alturapunto
            plt.plot((xigtw,xid),(yigtw,yid),
                     color = 'green',
                     linestyle='dotted')
            plt.scatter(xid ,yid)
            plt.annotate(unpunto,(xid,yid), rotation=45)
    plt.title('Linea :'+LineaObs)
    plt.xlabel('distancia')
    plt.ylabel('altura')
    plt.grid()

    # plt.axis('equal')   
    return(puntos)

def revisa_fresnel(analiza_punto, archivo_ruta,alturaGw,alturapunto,
                  plantas,guarda=False,nFres =1,freq=915,muestras=41):
    LineaObs = archivo_ruta.strip('.csv').split('_')[1]

    planta_desde  = plantas[0]
    planta_hasta  = plantas[1]
    planta_altura = plantas[2]

    # lambda Fresnel
    lamb = 300/freq #300e6(m/s)/(freq*1e6)
    puntos = pd.read_csv(archivo_ruta,index_col='punto')

    # busca gateway
    i_Gw    = 0
    i_punto = 0
    n_Gtw = ''
    i = 0
    for unpunto in puntos.index:
        if unpunto.startswith('Gw'):
            n_Gtw = unpunto
        if unpunto==analiza_punto:
            i_punto = i
        i = i+1

    # perfil   
    xi = np.array(puntos['dist_'+n_Gtw])
    yi = np.array(puntos['altitud'])
    etiqueta = puntos.index

    # plantaciones
    conplantas = (xi > planta_desde) & (xi < planta_hasta)
    plantacion = yi[conplantas]+planta_altura

    # añade alturas de antenas
    yi_ant = np.copy(yi) + alturapunto
    yi_ant[i_Gw] = yi[i_Gw] + alturaGw

    # Zona de Fresnel
    if i_punto>0:
        
        # linea directa
        dy = yi_ant[i_punto]-yi_ant[i_Gw]
        dx = xi[i_punto]-xi[i_Gw]
        m = dy/dx
        alpha = np.arctan(m)
        yi_f = lambda x: m*x + yi_ant[i_Gw]

        xi_D = np.linspace(xi[i_Gw],xi[i_punto],
                           muestras)
        yi_D = yi_f(xi_D)
        dist_D = np.sqrt(dx**2+dy**2)

        # zona Fresnel 1
        F1_up = np.zeros(len(xi_D))
        F1_down = np.zeros(len(xi_D))
        for i in range(0,len(xi_D),1):
            d1 = xi_D[i]/np.cos(alpha)
            d2 = dist_D-d1
            F1 = np.sqrt(np.abs(nFres*lamb*d1*d2)/(d1+d2))

            xiu = xi_D[i]-F1*np.sin(np.abs(alpha))
            d1_u = xiu
            d2_u = dist_D-d1_u

            xid = xi_D[i]+F1*np.sin(np.abs(alpha))
            d1_d = xid
            d2_d = dist_D-d1_d
            # Fresnel, formula
            Fup = np.sqrt(np.abs(nFres*lamb*d1_u*d2_u)/(d1_u+d2_u))
            Fdown = np.sqrt(np.abs(nFres*lamb*d1_d*d2_d)/(d1_d+d2_d))
            
            # Linea directa +/- Fresnel
            F1_up[i]   = yi_D[i] + Fup*np.cos(alpha)
            F1_down[i] = yi_D[i] - Fdown*np.cos(alpha)
    else:
        print('punto no encontrado',analiza_punto)

    # SALIDA
    # print('Radio Fresnel: ',np.max(F1u[i]))

    # GRAFICA
    plt.fill_between(xi[conplantas],yi[conplantas],
                     plantacion,color='lightgreen',
                     label = 'plantas')
    plt.plot(xi,yi,label='perfil', color = "brown")
    # antenas
    plt.scatter(xi,yi_ant,label=
                'Dispositivos',color='blue')
    if i_punto>0:
        # Fresnel
        plt.plot(xi_D,yi_D,
                 label='Directo',color='orange',
                 linestyle = 'dashed')
        plt.plot(xi_D,F1_down,
                 label='Fresnel1',color='orange')
        plt.plot(xi_D,F1_up,
                 color='orange')
        
    for i in range(0,len(xi),1):
        plt.annotate(etiqueta[i],
                     (xi[i],yi_ant[i]),rotation=45)
    plt.xlabel('distancia')
    plt.ylabel('altura')
    plt.title('Línea :'+LineaObs+', Enlace: '+n_Gtw+'-'+etiqueta[i_punto])
    plt.legend()
    plt.grid()
    # plt.axis('equal')
    # plt.show()
    return(puntos)

3.2 Ubica los puntos en gráfica 2D o 3D

Se grafica la ubicación de los puntos usando las coordenadas este, norte y altura registradas con el GPS. Las gráficas se pueden realizar en 2D y en 3D, lo que ofrece una vista de la ubicación general de los puntos respecto a los gateways.

La vista en 2d es semejante a la que se encuentra en GoogleEarth

Instrucciones en Python

# Datos desde GPS archivo.csv
# 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 = 'Coordenadas/coord_puntos.csv'
arch_rotacion = 'Coordenadas/rotando3D.gif'

# Parametros de gráfica
alturaantena = 8
tipograf = '3D' # '2D','3D'
etiquetar = ['Gw03','LOS08','CD01','CD11','CAG1','MA10']
grupo   = ['LOS','Ref','M','C', ]
colores = ['grey','magenta','green','brown',]

# PROCEDIMIENTO
# leer coordenadas
ubica = pd.read_csv(arch_gpsrsm, index_col='punto')
ubica = pd.DataFrame(ubica)
n = len(ubica)

# colores y alturas
ubica['color'] = 'Red'
for unpunto in ubica.index:
    # revisa grupo
    for i in range(0,len(grupo),1):
        if unpunto.startswith(grupo[i]):
            ubica.at[unpunto,'color'] = colores[i]
    if unpunto.startswith('Gw'):
        ubica.at[unpunto,'altitud'] = ubica.at[unpunto,'altitud']+alturaantena
print (ubica['color'])
    
# SALIDA
print(ubica.head())

# Grafica 2D --------------------------
if tipograf == '2D':
    figura, grafica = plt.subplots()
    for unpunto in ubica.index:
        unalon = ubica['utm_este'][unpunto]
        unalat = ubica['utm_norte'][unpunto]
        unaalt = ubica['altitud'][unpunto]
        uncolor = ubica['color'][unpunto] 
        condicion = not(unpunto.startswith('Ref'))
        if condicion:
            grafica.scatter(unalon,unalat,
                            color = uncolor,
                            label = unpunto)
            if unpunto in etiquetar:
                grafica.annotate(unpunto,(unalon,unalat),
                                 color = uncolor,
                                 rotation=45)
    plt.xlabel('UTM_este')
    plt.ylabel('UTM_norte')
    plt.title('Ubicacion UTM')
    plt.grid()
    plt.show()

# Grafica 3D --------------------------
if tipograf == '3D':
    figura = plt.figure()
    grafica = figura.add_subplot(projection='3d')
    # grafica = plt.figure().add_subplot(projection='3d')
    for unpunto in ubica.index:
        unalon  = ubica['utm_este'][unpunto]
        unalat  = ubica['utm_norte'][unpunto]
        unaalt  = ubica['altitud'][unpunto]
        uncolor = ubica['color'][unpunto] 
        condicion = not(unpunto.startswith('Ref'))
        if condicion:
            grafica.scatter(unalon,unalat,unaalt,
                            color = uncolor,
                            label = unpunto)
        if unpunto in etiquetar:
            grafica.text(unalon,unalat,unaalt,unpunto)
    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.1. Coordenadas – Perfil y zona Fresnel

Observar el perfil del terreno en una ruta, y estimar la primera zona de Fresnel permite revisar algunas situaciones donde se presente bloqueo, dispersión, difracción de la señal en el enlace.

Los datos del perfil se obtienen del proceso anterior desde el siguiente archivo:

Granja2021_LOS.csv

En el proceso se añade al perfil la altura de las antenas en los dispositivos, la altura de la vegetación promedio y se genera la gráfica mostrada.

punto,altitud,altitud_gps,dist_Gw03,latitud,longitud,utm_este,utm_norte,utm_zLetra,utm_zNum,vegetacion,alt_antena
Gw03,65.0,68.3,0.0,-2.140439491218423,-79.96247385948875,615377.0,9763377.0,M,17,65.0,73.0
LOS01,64.0,66.16,24.74,-2.140222431275481,-79.96252795441376,615371.0,9763401.0,M,17,64.0,65.0
LOS02,59.0,66.16,75.93,-2.139770213855424,-79.96262716484223,615360.0,9763451.0,M,17,59.0,60.0
LOS03,58.0,66.16,99.64,-2.139553123308579,-79.96263630222327,615359.0,9763475.0,M,17,59.5,59.0
...

El resultado se guarda en un archivo.csv y en la gráfica.png

Instrucciones en Python

# Procesa perfiles y Zona Fresnel
# girni-fiec-espol Revisión: 20220724
# http://blog.espol.edu.ec/girni/

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os

# INGRESO
archivocsv  = 'Granja2021_LOS.csv'
carp_coord  = 'coordenadas'
carpeta_rsm = 'resultado_LOS'
LineaObs  = 'LOS'
segmentos = ['Gw03','LOS']
# Fresnel, selecciona dispositivo
analiza_punto = 'LOS09'

# vegetacion
plantas_entre  = [0, 0] #[90,250]
plantas_altura = 1.6

# antenas
alturaGw    = 8
alturapunto = 1

# parametros Fresnel
nFres = 1
freq  = 915 # 915MHz

# Gráficas
muestras = 41 # resolución
guarda = True

# PROCEDIMIENTO
procesoCompleto = 0

# ruta de archivo
if len(carp_coord)>0:
    carp_coord = carp_coord+'/'
archivo_ruta = carp_coord + archivocsv
if os.path.exists(archivo_ruta):
    puntostodos  = pd.read_csv(archivo_ruta,
                               index_col='punto')
    procesoCompleto = 1

if procesoCompleto == 1:
    # busca gateway
    Gw_nombre = '';
    for unpunto in puntostodos.index:
        if unpunto.startswith('Gw'):
            Gw_nombre = unpunto
    if len(Gw_nombre)>0:
        procesoCompleto = 2

if procesoCompleto == 2:
    #ordena por distancia
    puntostodos = puntostodos.sort_values(by=['dist_'+Gw_nombre])
    puntostodos['enruta']=0
    for unpunto in puntostodos.index:
        for elemento in segmentos:
            if unpunto.startswith(elemento):
                puntostodos.loc[[unpunto],['enruta']] = 1

    puntos = puntostodos[puntostodos['enruta']==1].copy()
    procesoCompleto = 3

if procesoCompleto == 3:
    # busca gateway y ubicación de punto observado
    i_punto = -1 ; i_Gw=-1
    i=0
    for unpunto in puntos.index:
        if unpunto.startswith('Gw'):
            i_Gw = i
        if unpunto==analiza_punto:
            i_punto = i
        i = i+1
    if i_punto>=0 and i_Gw>=0:
        procesoCompleto = 4
    
if procesoCompleto ==4:
    # plantación sector
    puntos['vegetacion'] = puntos['altitud']
    for unpunto in puntos.index:
        dist_punto = puntos['dist_'+Gw_nombre][unpunto]
        envegetacion = (dist_punto>=np.min(plantas_entre))
        envegetacion = envegetacion and (dist_punto<=np.max(plantas_entre))
        if envegetacion:
            puntos.at[unpunto,'vegetacion'] = puntos['altitud'][unpunto] + plantas_altura

    # añade alturas de antenas
    puntos['alt_antena'] = puntos['altitud'] + alturapunto
    puntos.at[Gw_nombre,'alt_antena'] = puntos['altitud'][Gw_nombre] + alturaGw

    procesoCompleto = 5

if procesoCompleto == 5:
    # lambda Fresnel
    lamb = 300/freq #300e6(m/s)/(freq*1e6)

    # Zona de Fresnel
    xi = puntos['dist_'+Gw_nombre]
    yi_ant = puntos['alt_antena']
    if i_punto>0:
        
        # linea directa
        dy = yi_ant[i_punto]-yi_ant[i_Gw]
        dx = xi[i_punto]-xi[i_Gw]
        m = dy/dx
        alpha = np.arctan(m)
        yi_f = lambda x: m*x + yi_ant[i_Gw]

        xi_D = np.linspace(xi[i_Gw],xi[i_punto],
                           muestras)
        yi_D = yi_f(xi_D)
        dist_D = np.sqrt(dx**2+dy**2)

        # zona Fresnel 1
        F1_up = np.zeros(len(xi_D))
        F1_down = np.zeros(len(xi_D))
        for i in range(0,len(xi_D),1):
            d1 = xi_D[i]/np.cos(alpha)
            d2 = dist_D-d1
            F1 = np.sqrt(np.abs(nFres*lamb*d1*d2)/(d1+d2))

            xiu = xi_D[i]-F1*np.sin(np.abs(alpha))
            d1_u = xiu
            d2_u = dist_D-d1_u

            xid = xi_D[i]+F1*np.sin(np.abs(alpha))
            d1_d = xid
            d2_d = dist_D-d1_d
            # Fresnel, formula
            Fup = np.sqrt(np.abs(nFres*lamb*d1_u*d2_u)/(d1_u+d2_u))
            Fdown = np.sqrt(np.abs(nFres*lamb*d1_d*d2_d)/(d1_d+d2_d))
            
            # Linea directa +/- Fresnel
            F1_up[i]   = yi_D[i] + Fup*np.cos(alpha)
            F1_down[i] = yi_D[i] - Fdown*np.cos(alpha)    

    procesoCompleto = 6


# SALIDA
if procesoCompleto == 6:
    print(puntos)
    # guarda el reporte en csv
    unreporte = carpeta_rsm+'/perfil_'+LineaObs+'.csv'
    if os.path.exists(carpeta_rsm):
        puntos.to_csv(unreporte)
    else:
        print('no existe carpeta:',carpeta_rsm)
else:
    print('proceso completado hasta:', procesoCompleto)
    print('Revisar errores en proceso:', procesoCompleto+1)
    print('   Procesos:')
    print('1: abrir archivo de datos ',archivocsv )
    print('2: encontrar gateway con nombre Gw##')
    print('3: ordenar tabla por distancias al gateway')
    print('4: encontrar punto a analizar:',analiza_punto)
    print('5: añadir altura de antenas y plantas')
    print('6: calcular zona de fresnel')

if procesoCompleto == 6:
    # GRAFICA
    figura,(grafica1,grafica2) = plt.subplots(2,1)

    # vegetacion
    con_veget = np.array(puntos['vegetacion']-puntos['altitud'])
    xiv = np.array(puntos['dist_'+Gw_nombre])[con_veget>0]
    perfil = np.array(puntos['altitud'])[con_veget>0]
    vegetacion = np.array(puntos['vegetacion'])[con_veget>0]

    grafica1.fill_between(xiv,perfil,vegetacion,
                          color='lightgreen')
    # perfil
    grafica1.plot(puntos['dist_'+Gw_nombre],
                  puntos['altitud'],label='Perfil',
                  color='brown')
        
    # antenas
    grafica1.scatter(puntos['dist_'+Gw_nombre],
                     puntos['alt_antena'],
                     color='blue')
    # etiquetas y lineas enlace
    xigw = puntos['dist_'+Gw_nombre][Gw_nombre]
    yigw = puntos['alt_antena'][Gw_nombre]
    for unpunto in puntos.index:
        xip = puntos['dist_'+Gw_nombre][unpunto]
        yip = puntos['alt_antena'][unpunto]
        grafica1.annotate(unpunto,(xip,yip),
                          rotation=45)
        grafica1.plot((xigw,xip),(yigw,yip),
                      color = 'green',
                      linestyle='dotted')
    untitulo = 'Ruta: '+LineaObs
    untitulo = untitulo + ', Antena: Gw '+str(alturaGw)
    untitulo = untitulo + 'm Dispositivo '+str(alturapunto)
    untitulo = untitulo + 'm'
    grafica1.set_title(untitulo)
    grafica1.set_ylabel('altura')
    grafica1.legend()
    grafica1.grid(True,linestyle='dotted',
                  axis='x',which='both')
    # plt.axis('equal')

    grafica2.fill_between(xiv,perfil,vegetacion,
                          color='lightgreen')
    # perfil
    grafica2.plot(puntos['dist_'+Gw_nombre],
                  puntos['altitud'],
                  color='brown')
        
    # antenas
    grafica2.scatter(puntos['dist_'+Gw_nombre],
                     puntos['alt_antena'],
                     color='blue')
    # etiquetas y lineas enlace
    xigw = puntos['dist_'+Gw_nombre][Gw_nombre]
    yigw = puntos['alt_antena'][Gw_nombre]
    for unpunto in puntos.index:
        if unpunto in [Gw_nombre, analiza_punto]:
            xip = puntos['dist_'+Gw_nombre][unpunto]
            yip = puntos['alt_antena'][unpunto]
            grafica2.annotate(unpunto,(xip,yip),
                              rotation=45)

    if i_punto>0:
        # Fresnel
        grafica2.plot(xi_D,yi_D,
                      label='linea enlace',
                      color='orange',
                      linestyle = 'dashed')
        grafica2.plot(xi_D,F1_down,
                      label='Fresnel1',color='orange')
        grafica2.plot(xi_D,F1_up,
                      color='orange')

    grafica2.set_xlabel('distancia')
    grafica2.set_ylabel('altura')
    grafica2.legend()
    grafica2.grid(True,linestyle='dotted',
                  axis='x',which='both')
    # plt.axis('equal')

    unarchivo = carpeta_rsm+'/perfil_'+LineaObs+'.png'
    if (guarda==True) and os.path.exists(carpeta_rsm):
        plt.savefig(unarchivo)

    plt.show()

Referencia: https://es.wikipedia.org/wiki/Zona_de_Fresnel

3. Coordenadas – GPS a utm y distancias a gateway

Las coordenadas de cada punto donde se ubica un dispositivo se obtienen usando un GPS.

Cada nodo se identifica con tres letras y dos dígitos por ejemplo usando cultivo y numeración ascendente al alejarse del gateway. Por ejemplo, para el cultivo de Cacao , radial «A» desde el gateway y punto en ubicación 04  el identificador de ubicación será CA04.

Las coordenadas se registran obteniendo el archivo kml del equipo GPS y contienen los valores de Latitud y Longitud. Por facilidad de uso y ubicación las coordenadas se usan en formato utm (Universal Transverse Mercator).

El archivo KML contiene los datos de los puntos, básicamente en el siguiente formato:

<Placemark>
        <name>CA04
        <styleUrl>#msn_pink-diamond</styleUrl>
        <Point>
                <coordinates>-79.96257337394334,
                             -2.139534989113919,
                             65.76000000000001</coordinates>
        </Point>
</Placemark>

Se adjunta un ejemplo del archivo para pruebas.

Granja2021_LOS.kml

El objetivo es generar un archivo con distancias al gateway y alturas del perfil del terreno semejante al ejemplo mostrado:

       altitud  altitud_gps  dist_Gw03   latitud   longitud  utm_este  utm_norte utm_zLetra utm_zNum
punto                                                                                               
Gw03      65.0        68.30       0.00 -2.140439 -79.962474  615377.0  9763377.0          M       17
LOS01     64.0        66.16      24.74 -2.140222 -79.962528  615371.0  9763401.0          M       17
LOS02     59.0        66.16      75.93 -2.139770 -79.962627  615360.0  9763451.0          M       17
LOS03     58.0        66.16      99.64 -2.139553 -79.962636  615359.0  9763475.0          M       17
LOS04     58.0        66.06     126.44 -2.139309 -79.962645  615358.0  9763502.0          M       17
...

En cada punto, la altitud_gps se puede mejorar usando el perfil del terreno desde google earth. Las instrucciones se muestran más adelante en «Perfil de alturas de terreno».

Coordenadas en utm

Para la conversión a utm, se usa la librería en Python con el mismo nombre. 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 a UTM desde latitud y longitud y viceversa.

Las coordenadas en formato UTM para el ejemplo dado, se encuentran ubicadas en la zona «17 M». Revisar para cada caso la zona

Perfil de alturas de terreno

Los datos de alturas del gps usado presentan variaciones irregulares de altura, por lo que para mejores valores se usa el perfil desde google earth usando las marcas de las coordenadas del archivo.kml

Para visualizar el perfil se agrega una línea como una «ruta» entre el punto del Gateway (Gw3) y el último punto de observación. Los datos del perfil se pueden observar al seleccionar la ruta y usar la opción de «mostrar perfil de elevación». La opción se obtiene con «click derecho» del mouse. Se tiene adicionalmente como referencia, las distancias desde el gateway obtenidas en el proceso de conversión a utm. Los valores leidos de alturas para cada punto se guardan en un archivo .csv

El archivo obtenido en el ejemplo de alturas se muestra en:

Granja2021_alturas.csv


Instrucciones en Python

El algoritmo integra los resultados de la conversión de coordenadas a utm, incluye las distancias al gateway e integra las alturas obtenidas de forma externa para tener como resultado el archivo.csv.

Granja2021_LOS.csv

Una revisión para generar el archivo es evitar sobreescribir un archivo anterior. Por lo que si se requiere procesar nuevamente los datos se pide borrar el archivo anterior antes de ejecutar el algoritmo.

El algoritmo se divide en procesos identificados por un número ascendente. La variable ‘procesoCompletado‘ permite dar indicaciones de avances en el algoritmo en caso que se produzcan errores como haber escrito el nombre del archivo de datos con errores, que el nombre del gateway debe empezar con «Gw», etc.

# Procesa puntos de coordenadas, a utm
# girni-fiec-espol Revision:20220724
# http://blog.espol.edu.ec/girni/ 

import geopandas as gpd
import fiona
import utm
import numpy as np
import pandas as pd
import os

# INGRESO
archivokml  = 'Granja2021_LOS.kml' 
carp_coord  = 'Coordenadas'

# archivo de alturas tomadas de google
arch_altura = 'Granja2021_alturas.csv'

# zona UTM
zonaNum   = 17
zonaLetra = 'M'

precision = 2

# PROCEDIMIENTO
procesoCompleto = 0

# Lectura de archivo .KML
if len(carp_coord)>0:
    carp_coord = carp_coord+'/'
archivo_ruta = carp_coord+archivokml
if os.path.exists(archivo_ruta):
    gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'
    puntos = gpd.read_file(archivo_ruta, driver='KML')
    puntos = puntos.sort_values(by=['Name'])
    
    procesoCompleto = 1

if procesoCompleto == 1:
    # tabla en UTM
    tabla = pd.DataFrame()
    n_Gw = '';
    for i in puntos.index:
        # en latitud, longitud y altura
        punto_nombre = puntos.loc[i]['Name']
        longitud = puntos.loc[i]['geometry'].x
        latitud  = puntos.loc[i]['geometry'].y
        altitud  = puntos.loc[i]['geometry'].z
        altitud  = np.round(altitud,precision)
        #convierte lat_lon en UTM
        coord_utm = utm.from_latlon(latitud,longitud,
                                    zonaNum,zonaLetra)
        utm_este  = np.round(coord_utm[0],precision)
        utm_norte = np.round(coord_utm[1],precision)
        if punto_nombre.startswith('Gw'):
            n_Gw = punto_nombre
            x1 = utm_este
            y1 = utm_norte
                # añade coordenadas UTM a tabla
        tabla.at[i,'punto']      = punto_nombre
        tabla.at[i,'utm_este']   = utm_este
        tabla.at[i,'utm_norte']  = utm_norte
        tabla.at[i,'utm_zNum']   = str(zonaNum)
        tabla.at[i,'utm_zLetra'] = zonaLetra
        tabla.at[i,'altitud']    = 0
        tabla.at[i,'altitud_gps']= altitud
        tabla.at[i,'latitud']    = latitud
        tabla.at[i,'longitud']   = longitud
    procesoCompleto = 2

if procesoCompleto==2:
    # gateway encontrado
    if len(n_Gw)>0:
        procesoCompleto = 3

if procesoCompleto == 3:
    # distancias a Gateways
    for i in puntos.index:
        dist = 0.0
        x2 = tabla.at[i,'utm_este']
        y2 = tabla.at[i,'utm_norte']
        dist = np.sqrt((x2-x1)**2 + (y2-y1)**2)
        dist = np.round(dist,precision)
        tabla.at[i,'dist_'+n_Gw] = dist

    # indice por nombre de punto
    tabla = tabla.set_index('punto')
    # ordena por indice
    tabla = tabla.sort_index(axis=1)

    procesoCompleto = 4

if procesoCompleto == 4:
    # añade alturas desde archivo
    if len(arch_altura)>0:
        arch_rutaA = carp_coord + arch_altura
        if os.path.exists(arch_rutaA):
            tablaAltura = pd.read_csv(arch_rutaA)
            tablaAltura = tablaAltura.set_index('punto')
            for cadapunto in tabla.index:
                if cadapunto in tablaAltura.index:
                    encontrada = tablaAltura.loc[cadapunto]['altitud']
                    tabla.at[cadapunto,'altitud'] = encontrada

            procesoCompleto = 5
        
# Para escritura de archivo
archivocsv   = archivokml.split('.')[0] + '.csv'
archivo_ruta = carp_coord+archivocsv

# SALIDA
if procesoCompleto == 5:
    print(tabla)
    # Escribe archivo: tabla de resultados
    if os.path.exists(archivo_ruta):
        print('El archivo existe: ',archivo_ruta)
        print('Esta acción borraría el archivo último actualizado')
        print('Puede borrarlo manualmente antes de reintentar')
    else:
        tabla.to_csv(archivo_ruta)
else:
    print('proceso completado hasta:', procesoCompleto)
    print('Revisar errores en proceso:', procesoCompleto+1)
    print('   Procesos:')
    print('1: abrir archivo .KML',archivokml)
    print('2: crear tabla de datos')
    print('3: encontrar gateway:', n_Gw)
    print('4: calcula distancia cada punto a gateway')
    print('5: añadir alturas desde archivo: ', arch_rutaA)

Referencia: Bidirectional UTM-WGS84 converter for python. https://pypi.org/project/utm/0.4.2/

2.3. Descriptor estadístico de un punto o dispositivo

Los datos registrados para un punto, tabulados en el proceso anterior requieren un análisis básico para obtener descriptores estadísticos que los representen. A partir del archivo ejemplo mostrado, se pueden obtener los siguientes resultados:

data_LOS22.csv


Por ejemplo, para analizar el nivel de señal (rssi_up, rssi_down) se requiere de los parámetos de media, desviación estándar, medias móviles, función de probabilidad de masa (pmf), etc. Los parámetros se pueden resumir en tablas y gráficas semejantes a la mostrada.

descriptor estadístico para los datos de la gráfica

descriptor
               rssi_up  rssi_down
count        98.000000   98.00000
mean        -86.336735  -77.44898
std           3.782608    1.14083
min         -97.000000  -81.00000
25%         -89.000000  -78.00000
50%         -86.000000  -77.00000
75%         -84.000000  -77.00000
max         -78.000000  -75.00000
error_trama   0.000000    0.00000

Para el proceso se requiere el nombre del punto, la carpeta o directorio donde se encuentra, la medida a observar (ej: rssi).

Un parámetro como la media puede se insuficiente para observar el comportamiento, por lo que se añaden las medias móviles.

Para la media móvil se debe indicar cada cuántas muestras se obtendrá el promedio, por lo que al usar varias se darán varias de ellas mediante la variable movAvg_cual.

Errores en una trama

Se consideran errores de medida a los valores fuera de un intervalo expresado en «medida_normal». Al indicar si un dato tiene un error es posible determinar la tasa de error en de trama, para el canal de subida y de bajada.

Registros a procesar

Los registros que se procesarán serán libres de errores de trama, donde se calculan la media, desviación estándar, mínimo, máximo. Los resultados son los mostrados en la tabla.

Funciones de Probabilidad de masa (pmf)

Para cada valor de medida (rssi), se obtiene la frecuencia relativa que se presentan en una gráfica pmf.

Archivos de resultado

Los resultados se guardan en archivos.csv y graficos.pmf para el registro del análisis de cada punto. Loa archivos se agrupan en una carpeta de resultados para su posterior análisis en conjunto.

describe_LOS22.csv


Instrucciones en Python

# Descriptores estadisticos de los datos.csv
# de un dispositivo, revisa media, desviació estándar, pmf
# graba unreporte.csv con pandas y genera gráficas
# http://blog.espol.edu.ec/girni/

import numpy as np
import json
import pandas as pd
import datetime as dt
import os
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import matplotlib.units as munits

# INGRESO
cualpunto   = "data_m_LOS22.csv"
carpeta_rsm = "resultado_Test"

medida         = "rssi"
medida_unidad  = "dBm"
medida_normal  = [-250,-1]
medida_grafica = [-100,-60]

movAvg_cual  = [8,16] #cada cuantas muestras
movAvg_color = ['lightgreen','orange']
precision   = 2
guarda      = True

# PROCEDIMIENTO
unmodelo = cualpunto.strip('.csv').split('_')[1]
unubicado   = cualpunto.strip('.csv').split('_')[2]
codigopunto = unmodelo+'_'+unubicado

# Leer archivo
archivopunto = carpeta_rsm + '/' + cualpunto
tabla = pd.read_csv(archivopunto)
tabla = tabla.drop(columns='Unnamed: 0')
tabla = pd.DataFrame(tabla)

fechaformato = "%Y-%m-%d %H:%M:%S.%f"
# medida intervalo
medida_min = np.min(medida_normal)
medida_max = np.max(medida_normal)

# fechas series a datetime
tabla['publishedAt'] = pd.to_datetime(tabla['publishedAt'],
                                      format=fechaformato)
tabla['created'] = pd.to_datetime(tabla['publishedAt'],
                                  format=fechaformato)

# revisa errores de medida
tabla["error_up"]   = 0
tabla["error_down"] = 0
for undato in tabla['publishedAt'].keys():   
    medida_up = tabla[medida+'_up'][undato]
    enrango = (medida_up>=np.min(medida_normal))
    enrango = (enrango and medida_up<=np.max(medida_normal))
    if not(enrango):
        tabla.at[undato,"error_up"] = 1
    
    medida_down = tabla[medida+'_down'][undato]
    enrango = (medida_down>=np.min(medida_normal))
    enrango = (enrango and medida_down<=np.max(medida_normal))
    if not(enrango):
        tabla.at[undato,"error_down"] = 1      

# tasa error trama
leidos = len(tabla)
if leidos > 0:
    error_up   = np.sum(tabla['error_up'])
    error_up   = error_up/leidos
    error_down = np.sum(tabla['error_down'])
    error_down = error_down/leidos

# descriptor estadístico, datos sin errores
condicion_up = (tabla['error_up']==0)
condicion_down = (tabla['error_down']==0)

medida_up = tabla[condicion_up][medida+'_up']
describe_up = medida_up.describe()
describe_up['error_trama'] = error_up

medida_down = tabla[condicion_down][medida+'_down']
describe_down = medida_down.describe()
describe_down['error_trama'] = error_down

descriptor = describe_up.copy()
descriptor = pd.concat([descriptor,describe_down],axis=1)
descriptor['dispositivo'] = tabla['dispositivo'][0]

# función de probabilidad de masa pmf
def medida_pmf(valores,undescriptor):
    pmin   = np.min(valores)
    pmax   = np.max(valores)
    tramo  = int(pmax-pmin)
    conteo = np.zeros(tramo+1,dtype=int)
    intervalo = np.arange(pmin,pmax+1,1)
    for valor in valores:
        donde = np.where(intervalo == valor)
        conteo[donde] = conteo[donde] + 1
    freq_relativa = np.array(conteo)/np.sum(conteo)
    unpmf = {'intervalo' : list(intervalo),
             'freq_relativa' : list(freq_relativa)}
    return(unpmf)

pmf_up   = medida_pmf(medida_up,describe_up)
pmf_down = medida_pmf(medida_down,describe_down)

pmf_punto = {'pmf':{'pmf_up'   : pmf_up,
                    'pmf_down' : pmf_down}}
pmf_punto = pd.DataFrame(pmf_punto)
pmf_punto['dispositivo'] = tabla['dispositivo'][0]

# Para gráficas
# medias moviles en movAvg_cual[]
serie_up  = pd.Series(medida_up)
movAvg_up_mean = []
movAvg_up_std = []
m = len(movAvg_cual)
for j in range(0,m,1):
    k = movAvg_cual[j]
    movAvg_up_mean.append(list(serie_up.rolling(k).mean()))
    movAvg_up_std.append(list(serie_up.rolling(k).std()))
    
serie_down = pd.Series(medida_down)
movAvg_down_mean = []
movAvg_down_std = []
for j in range(0,m,1):
    k = movAvg_cual[j]
    movAvg_down_mean.append(list(serie_down.rolling(k).mean()))
    movAvg_down_std.append(list(serie_down.rolling(k).std()))

movAvgData ={'movAvg_cual'   : movAvg_cual,
             'movAvg_color'  : movAvg_color,
             'movAvg_up_mean'  : movAvg_up_mean,
             'movAvg_down_mean': movAvg_down_mean,
             'movAvg_up_std'   : movAvg_up_std,
             'movAvg_down_std' : movAvg_down_std
             }

grafData ={'codigopunto' : codigopunto,
           'medida' : medida,
           'precision': precision,
           'medida_unidad' : medida_unidad,
           'medida_grafica': medida_grafica
           }

# SALIDA --------------------------
print('descriptor')
print(descriptor)

# guarda el reporte en csv
unarchivo = carpeta_rsm+'/describe_'+codigopunto+'.csv'
descriptor.to_csv(unarchivo)

unarchivo = carpeta_rsm+'/pmf_'+codigopunto+'.json'
pmf_punto.to_json(unarchivo)

unarchivo = carpeta_rsm+'/movavg_'+codigopunto+'.json'
with open(unarchivo, "w") as outfile:
    json.dump(movAvgData, outfile)

unarchivo = carpeta_rsm+'/grfdata_'+codigopunto+'.json'
with open(unarchivo, "w") as outfile:
    json.dump(grafData, outfile)


# GRAFICA -----
def graf_puntos_serie(tabla,descriptor,movAvgData,grafData):
    ''' grafica la serie de tiempo de cada punto
        añade medias móviles
    '''
    # ajuste de formato de fecha para eje x
    converter = mdates.ConciseDateConverter()
    munits.registry[np.datetime64] = converter
    munits.registry[dt.date] = converter
    munits.registry[dt.datetime] = converter

    # datos para grafica
    precision   = grafData['precision']
    medida = grafData['medida']
    medida_unidad  = grafData['medida_unidad']
    medida_grafica = grafData['medida_grafica']
    
    movAvg_cual  = movAvgData['movAvg_cual']
    movAvg_color = movAvgData['movAvg_color']
    
    media_up    = descriptor[medida+'_up']['mean']
    std_up      = descriptor[medida+'_up']['std']
    media_down  = descriptor[medida+'_down']['mean']
    std_down    = descriptor[medida+'_down']['std']

    # ajuste de intervalo eje y
    y_min = np.min([np.min(medida_grafica),
                    media_up - 2*std_up,
                    media_down - 2*std_down])
    y_max = np.max([np.max(medida_grafica),
                    media_up + 2*std_up,
                    media_down + 2*std_down])
    
    # selecciona sin error
    condicion_up = (tabla['error_up']==0)
    condicion_down = (tabla['error_down']==0)

    # grafica
    fig_serie,(graf_up,graf_down) = plt.subplots(2,1)
    
    # medida_up -----
    graf_up.plot(tabla[condicion_up]['publishedAt'],
                 tabla[condicion_up][medida+'_up'],
                 color='blue',marker ='.',
                 linestyle='')
    
    # medida_up, medias y std
    etiq_up = str(np.round(media_up,precision))+' +/- '
    etiq_up = etiq_up + str(np.round(std_up,precision))
    graf_up.axhline(media_up,
                    color='blue',label=etiq_up)
    graf_up.axhline(media_up-std_up,
                    color='blue',linestyle='dotted')
    graf_up.axhline(media_up+std_up,
                    color='blue',linestyle='dotted')
    
    # medida_up, medias móviles
    m = len(movAvg_cual)
    for j in range(0,m,1):
        k = str(movAvg_cual[j])
        graf_up.plot(tabla[condicion_up]['publishedAt'],
                     movAvgData['movAvg_up_mean'][j],
                     label='movAvg_'+k,
                     color=movAvg_color[j])
    
    graf_up.set_ylim(y_min,y_max)
    graf_up.set_ylabel(medida+'_up ('+medida_unidad+')',
                       color='blue')
    graf_up.legend()
    graf_up.grid(True,linestyle='dotted',
                 axis='x',which='both')

    # medida_down -------
    graf_down.plot(tabla[condicion_down]['publishedAt'],
                   tabla[condicion_down][medida+'_down'],
                   color='brown',marker ='.',
                   linestyle='')

    # medida_down, medias y std
    etiq_down = str(np.round(media_down,precision))+' +/- '
    etiq_down = etiq_down + str(np.round(std_down,precision))
    graf_down.axhline(media_down,
                      color='brown',label=etiq_down)
    graf_down.axhline(media_down+std_down,
                      color='brown',linestyle='dotted')
    graf_down.axhline(media_down-std_down,
                      color='brown',linestyle='dotted')
    
    # medida_down, medias moviles
    for j in range(0,m,1):
        k = str(movAvg_cual[j])
        graf_down.plot(tabla[condicion_down]['publishedAt'],
                       movAvgData['movAvg_down_mean'][j],
                       label='movAvg_'+k,
                       color=movAvg_color[j])
    
    graf_down.set_ylim(y_min,y_max)
    graf_down.set_xlabel('fecha')
    graf_down.set_ylabel(medida+'_down ('+medida_unidad+')',
                         color='brown')
    graf_down.legend()
    graf_down.grid(True,linestyle='dotted',
                   axis='x', which='both')
    graf_up.set_title('Serie: '+grafData['codigopunto']+' '+ medida)
    plt.tight_layout()
    return(fig_serie)

def graf_puntos_pmf(pmf_punto,descriptor,grafData):
    ''' grafica función de probabilida de masa
        para cada punto, media +/- std
    '''
    # datos para grafica
    x_pmfup   = pmf_punto['pmf']['pmf_up']['intervalo']
    y_pmfup   = pmf_punto['pmf']['pmf_up']['freq_relativa']
    x_pmfdown = pmf_punto['pmf']['pmf_down']['intervalo']
    y_pmfdown = pmf_punto['pmf']['pmf_down']['freq_relativa']
    
    precision   = grafData['precision']
    medida = grafData['medida']
    medida_unidad  = grafData['medida_unidad']
    medida_grafica = grafData['medida_grafica']
    
    media_up   = descriptor[medida+'_up']['mean']
    std_up     = descriptor[medida+'_up']['std']
    media_down = descriptor[medida+'_down']['mean']
    std_down   = descriptor[medida+'_down']['std']

    prob_max = 0.40
    # ajuste de intervalo eje y
    y_min = np.min([np.min(medida_grafica),
                    media_up - 2*std_up,
                    media_down - 2*std_down])
    y_max = np.max([np.max(medida_grafica),
                    media_up + 2*std_up,
                    media_down + 2*std_down])
    # grafica
    fig_pmf,graf_pmf = plt.subplots()
    etiq_up = str(np.round(media_up,precision)) +' +/- '
    etiq_up = etiq_up + str(np.round(std_up,precision))
    graf_pmf.plot(x_pmfup,y_pmfup,
                  label='media_up '+etiq_up,
                  color='blue')
    graf_pmf.axvline(media_up,color='blue')
    graf_pmf.axvline(media_up+std_up,
                     linestyle='dotted',color='blue')
    graf_pmf.axvline(media_up-std_up,
                     linestyle='dotted',color='blue')

    etiq_down = str(np.round(media_down,precision))+' +/- '
    etiq_down = etiq_down + str(np.round(std_down,precision))
    graf_pmf.plot(x_pmfdown,y_pmfdown,
                  label='media_down '+etiq_down,
                  color='brown')
    graf_pmf.axvline(media_down,color='brown')
    graf_pmf.axvline(media_down+std_down,
                     linestyle='dotted',color='brown')
    graf_pmf.axvline(media_down-std_down,
                     linestyle='dotted',color='brown')

    graf_pmf.set_title('pmf: '+grafData['codigopunto']+' '+medida)
    graf_pmf.set_xlim(y_min,y_max)
    graf_pmf.set_ylim(0,prob_max)
    graf_pmf.set_xlabel(medida+' ('+medida_unidad+')')
    graf_pmf.set_ylabel('frecuencia relativa')
    graf_pmf.legend()
    graf_pmf.grid(True,linestyle='dotted',
                  axis='x', which='both')
    return(fig_pmf)

def graf_puntos_std(tabla,descriptor,movAvgData,grafData):
    ''' grafica serie de std usando medias moviles
        para cada punto, media_std
    '''
    # ajuste de formato de fecha para eje x
    converter = mdates.ConciseDateConverter()
    munits.registry[np.datetime64] = converter
    munits.registry[dt.date] = converter
    munits.registry[dt.datetime] = converter
    
    # datos para grafica
    precision = grafData['precision']
    medida    = grafData['medida']
    medida_unidad = grafData['medida_unidad']

    movAvg_cual = movAvgData['movAvg_cual']
    movAvg_color = movAvgData['movAvg_color']

    # selecciona sin error
    condicion_up   = (tabla['error_up']==0)
    condicion_down = (tabla['error_down']==0)
    
    # ajuste de intervalo eje y
    y_min = 0
    y_max = np.max([2, 2*descriptor[medida+'_up']['std'],
                    2*descriptor[medida+'_down']['std']])
    # grafica
    fig_std,(graf_stdUp,graf_stdDown) = plt.subplots(2,1)
    
    # std up
    std_up = np.round(descriptor[medida+'_up']['std'],precision)
    graf_stdUp.axhline(std_up,label='std '+str(std_up),
                       color='blue')
    m = len(movAvg_cual)
    for j in range(0,m,1):
        k = str(movAvg_cual[j])
        graf_stdUp.plot(tabla[condicion_up]['publishedAt'],
                        movAvgData['movAvg_up_std'][j],
                        label='movAvg_'+k,
                        color=movAvg_color[j])
    graf_stdUp.set_ylim(y_min,y_max)
    graf_stdUp.set_ylabel('std_up ('+medida_unidad+')',
                          color='blue')
    graf_stdUp.legend()
    graf_stdUp.grid(True,linestyle='dotted',
                    axis='x', which='both')
    graf_stdUp.set_title('std: '+grafData['codigopunto']+' '+ medida)

    # std down
    std_down = np.round(descriptor[medida+'_down']['std'],precision)
    graf_stdDown.axhline(std_down,label='std '+str(std_down),
                         color='brown')
    for j in range(0,m,1):
        k = str(movAvg_cual[j])
        graf_stdDown.plot(tabla[condicion_down]['publishedAt'],
                          movAvgData['movAvg_down_std'][j],
                          label='movAvg_'+k,color=movAvg_color[j])
    graf_stdDown.set_ylim(y_min,y_max)
    graf_stdDown.set_xlabel('fecha')
    graf_stdDown.set_ylabel('std_down ('+medida_unidad+')',
                            color='brown')
    graf_stdDown.legend()
    graf_stdDown.grid(True,linestyle='dotted',
                      axis='x', which='both')
    plt.tight_layout()
    return(fig_std)

fig_serie = graf_puntos_serie(tabla,descriptor,movAvgData,grafData)
fig_pmf   = graf_puntos_pmf(pmf_punto,descriptor,grafData)
fig_std   = graf_puntos_std(tabla,descriptor,movAvgData,grafData)

if guarda==True:
    unarchivo = carpeta_rsm+'/serie_'+codigopunto+'.png'
    fig_serie.savefig(unarchivo)
    unarchivo = carpeta_rsm+'/pmf_'+codigopunto+'.png'
    fig_pmf.savefig(unarchivo)
    unarchivo = carpeta_rsm+'/std_'+codigopunto+'.png'
    fig_std.savefig(unarchivo)
plt.show()

Referencia: pmf, cdf en una señal de sonido. http://blog.espol.edu.ec/estg1003/senal-de-sonido-pmf-cdf/

2.2. Lectura de registros en states.json de un Punto con Python

El procesamiento de los datos históricos de «un sensor» en una ubicación se realiza con Python desde el archivo generado en la sección anterior. Esta forma se obtiene un archivo nuevo y de menor tamaño para procesar con los datos de solo el dispositivo que se propone analizar.

states20211022.json

El formato de la tabla es el mismo que el observado en «DB Browser».

Los campos a usar son principalmente:

  • «entity_id» para la selección de un sensor,
  • «attributes» donde previamente se configuró para guardar  el mensaje mqtt de cada transmisión
  • «created» que contine la fecha y hora de creación del registro. referenciada a GMT=0.

Parametros guardados en attributes

El campo «attributes» de cada registro contiene varios parámetros posibles a usar en el formato del mensaje. Un ejemplo del contenido se indica en el recuadro en rojo mostrado en lala figura anterior, y se presenta a continuación.

{"applicationID":"4",
 "applicationName":"Leer_RssiSNR",
 "deviceName":"cc11Pruebas",
 "devEUI":"pT7GFa7ePxE=",
 "rxInfo":[{"gatewayID":"uCfr//4dylc=",
            "time":null,"timeSinceGPSEpoch":null,
            "rssi":-63,
            "loRaSNR":10.8,
            "channel":3,
            "rfChain":0,
            "board":0,
            "antenna":0,
            "location":{"latitude":-2.1418242484279535,
                        "longitude":-79.96711134910585,
                        "altitude":20,
                        "source":"UNKNOWN",
                        "accuracy":0},
            "fineTimestampType":"NONE",
            "context":"CmMVZA==",
            "uplinkID":"I7NWhENCQyaFtYoZIDw3KQ==",
            "crcStatus":"CRC_OK"}],
 "txInfo":{"frequency":902900000,
           "modulation":"LORA",
           "loRaModulationInfo":{"bandwidth":125,
                                 "spreadingFactor":10,
                                 "codeRate":"4/5",
                                 "polarizationInversion":false}
           },
 "adr":true,
 "dr":0,
 "fCnt":778,
 "fPort":4,
 "data":"KwcK9A8=",
 "objectJSON":"{\"Down_datarate\":10,\"Down_rssi\":-43,\"Down_snr\":7,\"bateria_V\":4.08}",
 "tags":{},
 "confirmedUplink":true,
 "devAddr":"AVqCdA==",
 "publishedAt":"2021-10-21T16:28:36.265518462Z",
 "unit_of_measurement":"dBm",
 "friendly_name":"rssi_up_cc11"
 }

El formato del ejemplo mostrado se realizó con el editor de Python para facilidad de lectura,

Procesamiento de datos para un sensor

El procesamiento de cada registro se realiza separando los campos a usar, en una tabla los campos representan las  columnas a usar.

La selección del dispositivos se realiza comparando con «entity_id«.

El intervalo de fechas permite usar ventanas de tiempo de observación. La fecha y hora registrada se encuentra referenciadas a GMT=0, dado que en el pais donde se realizan las muestras usa GMT=-5, se realiza un desplazamiento en las horas para facilitar la identificación de las franjas horarias.

Los identificadores de gateways se simplifican con el diccionario gateway, que deben ser previamente identificados buscando en el archivo.json como «gatewayID». También se revisa el número de trama (fCnt) para registrar las recepción de tramas repetidas en cada gateway, registrando los valores en una lista, para posteriormente, seleccionar el valor mayor, menor, promedio).

Los parámetros de recepción en el gateway se encuentran en «rx_datos» en el siguiente orden: rssi, loRaSNR, channel, rfChain, frequency

El resultado de la selección de los datos se guarda en un nuevo archivo en formato.csv para facilitar la lectura de los datos con otros programas.

Otra forma de guardar los datos en en formato.json, que puede ser leida con la librería Pandas.

Instrucciones en Python para tabla CSV

# Lectura de archivo.json hacia una tabla0
# selecciona registros de "un sensor"
# graba unreporte.csv con pandas
# http://blog.espol.edu.ec/girni/

import numpy as np
import json
import pandas as pd
import datetime as dt
import os

# INGRESO
# archivo de entrada
unarchivoDB = "states20211022.json"
carpeta_DB  = "BaseDatosJSON"

# fecha intervalo
fechainicio = "2021-10-21 11:00:00.0"
fechafin    = "2021-10-22 11:50:00.0"
zonaGMT     = -5

unsensor    = "sensor.rssi_up_cc12"
modelo_disp = "modulo" # capsula, modulo desarrollo
ubicado     = "LOS22"
carpeta_rsm = "resultado_Test"

# gatewayID simplifica identificador 
gatewayDB   = {'uCfr//4dylc=':'Gw03'}

# PROCEDIMIENTO -----
campos = ['domain','entity_id','state','attributes','created']

# Lectura de archivo json, cambia a formato DataFrame
tabla0 = pd.DataFrame()
archivo_ruta  = carpeta_DB + '/' + unarchivoDB
archivoexiste = os.path.exists(archivo_ruta)
if archivoexiste:
    tabla0 = pd.read_json(archivo_ruta)
    tabla0 = pd.DataFrame(tabla0,columns=campos)
else:
    print(' ERROR: NO se encuentra el archivo...')
    print('        revise el nombre de archivo. ')

# Revisa cada registro 
leidos = 0
tabla  = {}
if archivoexiste:
    # Intervalo de fecha como datetime
    hora_desplaza  = dt.timedelta(hours = abs(zonaGMT))
    fechaformatoDB = '%Y-%m-%dT%H:%M:%S.%f'
    fechaformato   = '%Y-%m-%d %H:%M:%S.%f'
    fechatxt     = fechainicio[0:np.min([26,len(fechainicio)])]
    fechainicio  = dt.datetime.strptime(fechatxt,fechaformato)
    fechatxt     = fechafin[0:np.min([26,len(fechafin)])]
    fechafin     = dt.datetime.strptime(fechatxt,fechaformato)

    # datos de trama
    for cadaregistro in tabla0.index:
        unatrama   = tabla0['attributes'][cadaregistro]
        trama_mqtt = json.loads(unatrama)
        
        # selecciona registro por sensor, en fecha y con datos
        cualsensor = (tabla0['entity_id'][cadaregistro] == unsensor)
        
        unafecha = tabla0['created'][cadaregistro]
        unafecha = unafecha[0:np.min([26,len(unafecha)])]
        unafecha = dt.datetime.strptime(unafecha,fechaformato)
        unafecha = unafecha - hora_desplaza
        
        enfecha  = (unafecha>=fechainicio) and (unafecha<=fechafin)
        condatos = 'applicationID' in trama_mqtt.keys()
        
        selecciona = cualsensor and condatos and enfecha
        if (selecciona == True):        
            # datos en la mensaje MQTT
            publishedAt = trama_mqtt["publishedAt"]
            publishedAt = publishedAt[0:np.min([26,len(publishedAt)])]
            publishedAt = dt.datetime.strptime(publishedAt,fechaformatoDB)
            publishedAt = publishedAt - hora_desplaza
            
            friendly_name = trama_mqtt["friendly_name"]
            deviceName    = trama_mqtt["deviceName"]
            dispositivo   = friendly_name.split('_')[2]
            
            dr      = trama_mqtt["dr"]
            fCnt    = trama_mqtt["fCnt"]
            fPort   = trama_mqtt["fPort"]
            freq_tx = trama_mqtt["txInfo"]["frequency"]
            bandwidth = trama_mqtt["txInfo"]["loRaModulationInfo"]["bandwidth"]
            spreadingFactor = trama_mqtt["txInfo"]["loRaModulationInfo"]["spreadingFactor"]

            objectJSON = trama_mqtt["objectJSON"]
            datosensor = json.loads(objectJSON)
                 
            datostrama = {"publishedAt": publishedAt,
                          "dispositivo": dispositivo,
                          "fCnt": fCnt}
            
            # revisa gateway en la red LoRaWan
            tamano = len(trama_mqtt["rxInfo"])
            i = 0
            while i<tamano:
                gatewayID = trama_mqtt["rxInfo"][i]["gatewayID"]
                rssi      = trama_mqtt["rxInfo"][i]["rssi"]
                loRaSNR   = trama_mqtt["rxInfo"][i]["loRaSNR"]
                channel   = trama_mqtt["rxInfo"][i]["channel"]
                rfChain   = trama_mqtt["rxInfo"][i]["rfChain"]
                gtwNum    = gatewayDB[gatewayID]
                
                # registra en tabla, incluyendo tramas repetidas
                datostrama['rssi_up']    = rssi
                datostrama['snr_up']     = loRaSNR
                datostrama['channel_up'] = channel
                datostrama['rfChain_up'] = rfChain
                datostrama['gtw_rx']     = gtwNum
                i = i + 1

            # dato del sensor
            equivale = {'Down_datarate': 'dr_down',
                        'Down_rssi'    : 'rssi_down',
                        'Down_snr'     : 'snr_down',
                        'bateria_V'    : 'bateria_V'}
            for undato in datosensor:
                if undato in equivale.keys():
                    datostrama[equivale[undato]] = datosensor[undato]
                else:
                    datostrama[undato] = datosensor[undato]

            # datos, otros valores
            datostrama["frequency"] = freq_tx
            datostrama["bandwidth"] = bandwidth
            datostrama["spreadingFactor"] = spreadingFactor
            datostrama["fPort"]   = fPort
            datostrama["dr"]      = dr
            datostrama["created"] = unafecha

            leidos = leidos + 1
            
            # revisa registro repetido
            repetido = False
            if leidos>1:
                trama_num  = (fCnt == tabla[leidos-1]['fCnt'])
                fecha_pub  = (publishedAt == tabla[leidos-1]['publishedAt'])
                gtwNum_rep = (gtwNum == tabla[leidos-1]['gtw_rx'])
                repetido   = (trama_num and fecha_pub and gtwNum_rep)
            if not(repetido):
                tabla[leidos] = datostrama
            else:
                leidos = leidos - 1

    # convierte diccionario a DataFrame
    if len(tabla)>0:
        tabla = pd.DataFrame(tabla)
        tabla = tabla.T
        # añade ubicación y modelo (capsula o modulo desarrollo)
        tabla.insert(1,'ubicado',ubicado)
        tabla.insert(2,'modelo_disp',modelo_disp)
        
# directorio de resultados
if len(carpeta_rsm) > 0:
    if not(os.path.exists(carpeta_rsm)):
        os.makedirs(carpeta_rsm)
    carpeta_rsm = carpeta_rsm + '/'

# SALIDA
print('registros leidos:     ', len(tabla0))
print('registros procesados: ', leidos)
print(tabla)

# guarda el reporte en csv
if len(tabla)>0:
    unreporte = carpeta_rsm+'data_'+ubicado+'.csv'
    tabla.to_csv(unreporte)

el resultado de puede visualiza como:

registros leidos:      7090
registros procesados:  98
                  publishedAt ubicado modelo_disp dispositivo  ... spreadingFactor fPort dr                    created
1  2021-10-21 11:34:07.079188   LOS22      modulo        cc12  ...              10     4  0 2021-10-21 11:34:07.105036
2  2021-10-21 11:49:07.862255   LOS22      modulo        cc12  ...              10     4  0 2021-10-21 11:49:07.876219
3  2021-10-21 12:04:07.931128   LOS22      modulo        cc12  ...              10     4  0 2021-10-21 12:04:07.943781
4  2021-10-21 12:19:08.489097   LOS22      modulo        cc12  ...              10     4  0 2021-10-21 12:19:08.504848
5  2021-10-21 12:34:08.479274   LOS22      modulo        cc12  ...              10     4  0 2021-10-21 12:34:08.494229
..                        ...     ...         ...         ...  ...             ...   ... ..                        ...
94 2021-10-22 10:49:46.742877   LOS22      modulo        cc12  ...              10     4  0 2021-10-22 10:49:46.756601
95 2021-10-22 11:04:47.636628   LOS22      modulo        cc12  ...              10     4  0 2021-10-22 11:04:47.652299
96 2021-10-22 11:19:48.310535   LOS22      modulo        cc12  ...              10     4  0 2021-10-22 11:19:48.331143
97 2021-10-22 11:34:49.209626   LOS22      modulo        cc12  ...              10     4  0 2021-10-22 11:34:49.224717
98 2021-10-22 11:49:49.638194   LOS22      modulo        cc12  ...              10     4  0 2021-10-22 11:49:49.651398

[98 rows x 20 columns]
>>>

y el archivo resultante se presenta como:

data_LOS22.csv


Instrucciones en Python para tabla JSON

# Lectura de archivo.json hacia una tabla1
# selecciona registros de "un sensor"
# graba unreporte.json con pandas

import numpy as np
import json
import pandas as pd
import datetime as dt

# INGRESO

# archivo de entrada
unarchivo   = 'states20211022.json'

fechainicio = '2021-10-21 19:00:00.0'
fechafin    = '2021-10-22 17:10:00.0'

unsensor    = 'sensor.rssi_up_cc22'
# archivo de salida
unreporte   = 'tramas20211022_cc12.json'

# PROCEDIMIENTO -----
# columnas o campos a usar
campos    = ['domain','entity_id','state',
             'attributes','created']
# gatewayID simplifica identificador 
gateway   = {'uCfr//5zvhk=':'Gw01',
             'uCfr//4dylc=':'Gw03'}

# Lectura de archivo json
tabla0 = pd.read_json(unarchivo)
# Cambia a formato DataFrame
tabla0 = pd.DataFrame(tabla0,columns=campos)

# Revisa cada registro 
leidos = 0
tabla  = {}

# Intervalo de fecha como datetime
formato     = '%Y-%m-%d %H:%M:%S.%f'
fechainicio = fechainicio[0:np.min([26,len(fechainicio)])]
fechainicio = dt.datetime.strptime(fechainicio,formato)
fechafin    = fechafin[0:np.min([26,len(fechafin)])]
fechafin    = dt.datetime.strptime(fechafin,formato)

# datos de trama
for cadaregistro in tabla0.index:
    unatrama   = tabla0['attributes'][cadaregistro]
    trama_mqtt = json.loads(unatrama)
    
    # selecciona registro por sensor, en fecha y con datos
    cualsensor = (tabla0['entity_id'][cadaregistro] == unsensor)
    
    unafecha = tabla0['created'][cadaregistro]
    unafecha = unafecha[0:np.min([26,len(unafecha)])]
    unafecha = dt.datetime.strptime(unafecha,formato)
    
    enfecha  = (unafecha>=fechainicio) and (unafecha<=fechafin)
    condatos   = 'applicationID' in trama_mqtt.keys()
    
    selecciona = cualsensor and condatos and enfecha
    if (selecciona == True):
        leidos =  leidos + 1
        
        # datos en la mensaje MQTT
        publishedAt = trama_mqtt["publishedAt"]
        friendly_name = trama_mqtt["friendly_name"]
        deviceName    = trama_mqtt["deviceName"]
        dr    = trama_mqtt["dr"]
        fCnt  = trama_mqtt["fCnt"]
        fPort = trama_mqtt["fPort"]

        objectJSON  = trama_mqtt["objectJSON"]
        datosensor  = json.loads(objectJSON)
        
        frequency   = trama_mqtt["txInfo"]["frequency"]
        bandwidth   = trama_mqtt["txInfo"]["loRaModulationInfo"]["bandwidth"]
        spreadingFactor = trama_mqtt["txInfo"]["loRaModulationInfo"]["spreadingFactor"]
             
        datostrama = {"publishedAt":publishedAt,
                      "fPort":fPort, "dr": dr,"fCnt": fCnt,
                      "frequency":  frequency,
                      "bandwidth":  bandwidth,
                      "spreadingFactor": spreadingFactor,
                      "datosensor":datosensor}
        
        # revisa por cada gateway en la red LoRaWan
        tamano = len(trama_mqtt["rxInfo"])
        i = 0
        while i<tamano:
            gatewayID = trama_mqtt["rxInfo"][i]["gatewayID"]
            rssi      = trama_mqtt["rxInfo"][i]["rssi"]
            loRaSNR   = trama_mqtt["rxInfo"][i]["loRaSNR"]
            channel   = trama_mqtt["rxInfo"][i]["channel"]
            rfChain   = trama_mqtt["rxInfo"][i]["rfChain"]
            gtwNum    = gateway[gatewayID]
            i = i+1
            
            # registra en tabla, incluyendo tramas repetidas
            rx_datos = [rssi, loRaSNR,
                        channel, rfChain,frequency/1e6]  
            if not(gtwNum in datostrama.keys()):
                datostrama[gtwNum] = [rx_datos]
            else:
                if not(rx_datos in datostrama[gtwNum]):
                    datostrama[gwNum].append(rx_datos)
        dispositivo = friendly_name.split('_')[2]

        # revisa registro repetido
        disp_existe = dispositivo in tabla
        if not(disp_existe):
            tabla[dispositivo] = {leidos:datostrama}
        else:
            trama_num = (fCnt == tabla[dispositivo][leidos-1]['fCnt'])
            fecha_pub = (publishedAt == tabla[dispositivo][leidos-1]['publishedAt'])
            repetido  = (trama_num and fecha_pub)
            if not(repetido):
                tabla[dispositivo][leidos] = datostrama
            else:
                leidos = leidos - 1

# convierte diccionario en DataFrame
tabla = pd.DataFrame(tabla)

# SALIDA
print('registros leidos:     ',len(tabla0))
print('registros procesados: ', leidos)

# guarda el reporte en json
tabla.to_json(unreporte)

El procesamiento de los datos se puede realizar con librerias como Pandas.

En el procesamiento de registros, se descarta el registro de «join».

Una muestra del reporte es:

registros leidos:      7090
registros procesados:  98
>>> tabla.keys()
Index(['cc12'], dtype='object')
>>> tabla['cc12'].keys()
Int64Index([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
            18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
            35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
            52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68,
            69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85,
            86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98],
           dtype='int64')
>>> tabla['cc12'][2].keys()
dict_keys(['publishedAt', 'fPort', 'dr', 'fCnt',
 'frequency', 'bandwidth', 'spreadingFactor', 'datosensor',
 'Gw03'])
>>> 

>>> tabla['cc12'][2]['datosensor'].keys()
dict_keys(['Down_datarate', 'Down_rssi', 'Down_snr',
 'bateria_V'])

>>> tabla['cc12'][2]['Gw03']
[[-84, 8.5, 1, 0, 902.5]]

Esta forma se obtiene un nuevo y mas pequeño archivo para procesar con los datos de solo el dispositivo que se propone analizar.

tramas20211022_cc12.json

Referencia: Archivos.json – Pandas en blog de Fundamentos de Programación