Autor: Edison Del Rosario

3. LoRa Multipunto - HELTEC CubeCell SemiDuplex.ino

Como cada dispositivo contiene solo un radio, la comunicación en dos sentidos puede habilitarse en modo SemiDuplex. Las funciones de transmisión y recepción se alternan.

Resultados

TX Paquete "213110" , tamano 6 bytes , proximo en 8118 ms 
    RX Paquete "312139" , tamanio 6, Rssi -22 
TX Paquete "213111" , tamano 6 bytes , proximo en 8429 ms 
    RX Paquete "312140" , tamanio 6, Rssi -23 
TX Paquete "213112" , tamano 6 bytes , proximo en 9397 ms

Instrucciones

/* LoRa TRANSMITE/RECIBE Semi Duplex / Half-Duplex
 * Referencia: https://github.com/HelTecAutomation/ASR650x-Arduino
*/
#include "LoRaWan_APP.h"
#include "Arduino.h"

#ifndef LoraWan_RGB    // LED placa
#define LoraWan_RGB 0
#endif

// LoRa Parametros 
#define RF_FREQUENCY       915E6     // Hz
#define TX_OUTPUT_POWER    14        // dBm
#define LORA_BANDWIDTH     0         // [0: 125 kHz, 1: 250 kHz,
                                     //  2: 500 kHz, 3: Reserved]
#define LORA_SPREADING_FACTOR   7    // [SF7..SF12]
#define LORA_CODINGRATE         1    // [1: 4/5,  2: 4/6,
                                     //  3: 4/7,  4: 4/8]
#define LORA_PREAMBLE_LENGTH    8    // Same for Tx and Rx
#define LORA_SYMBOL_TIMEOUT     0    // Symbols
#define LORA_FIX_LENGTH_PAYLOAD_ON  false
#define LORA_IQ_INVERSION_ON        false
#define RX_TIMEOUT_VALUE            1000
#define BUFFER_SIZE                 30 // Define the payload size here
char txpacket[BUFFER_SIZE];          // cadena de caracteres
char rxpacket[BUFFER_SIZE];
static RadioEvents_t RadioEvents;

void OnTxDone( void );     // Tx completada
void OnTxTimeout( void );  // Tx fuera de tiempo
void OnRxDone( uint8_t *payload, uint16_t size, int16_t rssi, int8_t snr );

byte modoOperacion = 0; // 0:RX 1:TX
bool sleepMode = false;
int16_t txNumero;
int16_t Rssi,rxSize;

// Direcciones por dispositivo
byte dir_local   = 21; // Dispositivo envia
byte dir_destino = 31; // Dispositivo recibe

// tiempo entre Tx de datos o lecturas de sensor
long tiempo_antes     = 0;
long tiempo_intervalo = 7000;
long tiempo_espera = tiempo_intervalo + random(3000);

void setup() {
  Serial.begin(115200);
  txNumero=10; Rssi=0;

  RadioEvents.TxDone = OnTxDone;
  RadioEvents.TxTimeout = OnTxTimeout;
  RadioEvents.RxDone = OnRxDone;
  Radio.Init( &RadioEvents );

  Radio.SetChannel( RF_FREQUENCY );
  Radio.SetTxConfig( MODEM_LORA, TX_OUTPUT_POWER, 0, LORA_BANDWIDTH,
                     LORA_SPREADING_FACTOR, LORA_CODINGRATE,
                     LORA_PREAMBLE_LENGTH, LORA_FIX_LENGTH_PAYLOAD_ON,
                     true, 0, 0, LORA_IQ_INVERSION_ON, 3000 );
  Radio.SetRxConfig( MODEM_LORA, LORA_BANDWIDTH, LORA_SPREADING_FACTOR,
                     LORA_CODINGRATE, 0, LORA_PREAMBLE_LENGTH,
                     LORA_SYMBOL_TIMEOUT, LORA_FIX_LENGTH_PAYLOAD_ON,
                     0, true, 0, 0, LORA_IQ_INVERSION_ON, true );
  modoOperacion=0;  // 0:RX 1:TX 2:LOWPOWER
}

void loop(){
  // Intervalos entre mensajes
  long tiempo_ahora   = millis();
  long t_transcurrido = tiempo_ahora - tiempo_antes;
  if (t_transcurrido >= tiempo_espera){
    tiempo_antes = millis();
    tiempo_espera = tiempo_intervalo + random(3000);
    modoOperacion=1; // 0:RX 1:TX
  }else{
    modoOperacion=0; // 0:RX 1:TX
  }
  
  switch(modoOperacion) { // 0:RX 1:TX 2:LOWPOWER
    case 0:
      Radio.Rx( 0 );   
      modoOperacion=2;  //LOWPOWER;
      break;
    case 1:
      delay(500);
      transmiteMsg();
      modoOperacion=0;
      break;
    case 2:
      lowPowerHandler();
      modoOperacion=0;
      break;
    default:
      break;
  }
  turnOnRGB(0,0); // LED apaga
  delay(100);
  Radio.IrqProcess( );
}

Procedimiento de transmisión de paquete LoRa

void transmiteMsg( void ) {
  turnOnRGB(COLOR_SEND,0); // LED de placa
  // Paquete a transmitir
  sprintf(txpacket,"%d",dir_local);
  sprintf(txpacket+strlen(txpacket),"%d",dir_destino);
  sprintf(txpacket+strlen(txpacket),"%d",txNumero); //añade número paquete a txpacket

  // Mensaje a pantalla
  Serial.printf("\rTX Paquete \"%s\" , tamano %d bytes\r",
                txpacket, strlen(txpacket));
  Serial.printf("\r , proximo en %d ms \r\n",
                tiempo_espera);
  
  // Transmite paquete LoRa
  Radio.Send( (uint8_t *)txpacket, strlen(txpacket) );

  txNumero = txNumero + 1;  // cuenta paquete
  if (txNumero>=99){
    txNumero = 0;           // reinicia contador
  }
  modoOperacion=2;//LOWPOWER;
}

void OnTxDone( void ) {
  turnOnRGB(0,0);
  modoOperacion=0;
}

void OnTxTimeout( void ) {
  Radio.Sleep( );
  Serial.println("TX Timeout......");
  modoOperacion=1;
}

Procedimiento de recepción de paquete LoRa

void OnRxDone( uint8_t *payload, uint16_t size, int16_t rssi, int8_t snr ) {
  turnOnRGB(COLOR_RECEIVED,0);
  Rssi = rssi;     // nivel de recepcion
  rxSize = size;
  memcpy(rxpacket, payload, size );
  rxpacket[size]='\0';  //añade fin de cadena

  Radio.Sleep( );

  Serial.printf("\r    RX Paquete \"%s\" , tamanio %d, Rssi %d \r\n",
                rxpacket,rxSize,Rssi);
}

2022-12-08

2.2 LoRa Multipunto - HELTEC CubeCell Transmisor.ino

Resultados

TX Paquete "213110" , tamano 6 bytes , proximo en 8118 ms 
TX Paquete "213111" , tamano 6 bytes , proximo en 8429 ms 
TX Paquete "213112" , tamano 6 bytes , proximo en 9397 ms

Instrucciones

/* TRANSMITE de Mensajes Heltec Automation send communication test example
 * Referencia: https://github.com/HelTecAutomation/ASR650x-Arduino 
*/
#include "LoRaWan_APP.h"
#include "Arduino.h"

#ifndef LoraWan_RGB    // LED placa
#define LoraWan_RGB 0
#endif

// LoRa Parametros
#define RF_FREQUENCY       915E6     // Hz
#define TX_OUTPUT_POWER    14        // dBm
#define LORA_BANDWIDTH     0         // [0: 125 kHz, 1: 250 kHz,
                                     //  2: 500 kHz, 3: Reserved]
#define LORA_SPREADING_FACTOR   7    // [SF7..SF12]
#define LORA_CODINGRATE         1    // [1: 4/5,  2: 4/6,
                                     //  3: 4/7,  4: 4/8]
#define LORA_PREAMBLE_LENGTH    8    // Same for Tx and Rx
#define LORA_SYMBOL_TIMEOUT     0    // Symbols
#define LORA_FIX_LENGTH_PAYLOAD_ON  false
#define LORA_IQ_INVERSION_ON        false
#define RX_TIMEOUT_VALUE            1000
#define BUFFER_SIZE                 30 // Tamaño de paquete
char txpacket[BUFFER_SIZE];       // cadena de caracteres
char rxpacket[BUFFER_SIZE];
static RadioEvents_t RadioEvents;

byte  txNumero;
int16_t Rssi,rxSize;

// Direcciones por dispositivo
byte dir_local   = 31; // Dispositivo envia
byte dir_destino = 21; // Dispositivo recibe

// tiempo entre Tx de datos o lecturas de sensor
long tiempo_antes     = 0;
long tiempo_intervalo = 7000;
long tiempo_espera = tiempo_intervalo + random(3000);

void setup() {
  Serial.begin(115200);
  txNumero=10;  Rssi=0;

  Radio.Init( &RadioEvents );
  Radio.SetChannel( RF_FREQUENCY );
  Radio.SetTxConfig( MODEM_LORA, TX_OUTPUT_POWER, 0, LORA_BANDWIDTH,
                    LORA_SPREADING_FACTOR, LORA_CODINGRATE,
                    LORA_PREAMBLE_LENGTH, LORA_FIX_LENGTH_PAYLOAD_ON,
                    true, 0, 0, LORA_IQ_INVERSION_ON, 3000 ); 
}
void loop(){
  // Intervalos entre mensajes
  long tiempo_ahora   = millis();
  long t_transcurrido = tiempo_ahora - tiempo_antes;

  if (t_transcurrido >= tiempo_espera){
    transmiteMsg();  
    tiempo_antes = millis(); //actualiza tiempos
    tiempo_espera = tiempo_intervalo + random(3000);
  }
}

Procedimiento de transmisión de paquete LoRa

void transmiteMsg( void ) {
  turnOnRGB(COLOR_SEND,0); // LED de placa

  // Paquete a transmitir
  sprintf(txpacket,"%d",dir_local);
  sprintf(txpacket+strlen(txpacket),"%d",dir_destino);
  sprintf(txpacket+strlen(txpacket),"%d",txNumero); //añade número paquete a txpacket
  
  // Mensaje a pantalla
  Serial.printf("\rTX Paquete \"%s\" , tamano %d bytes\r",
                txpacket, strlen(txpacket));
  Serial.printf("\r , proximo en %d ms \r\n",
                tiempo_espera);

  // Transmite paquete LoRa
  Radio.Send( (uint8_t *)txpacket, strlen(txpacket) );
  
  txNumero = txNumero + 1;  // cuenta paquete
  if (txNumero>=99){
    txNumero = 0;           // reinicia contador
  }
}

void OnTxDone( void ) {
  turnOnRGB(0,0);
}

void OnTxTimeout( void ) {
  Radio.Sleep( );
  Serial.println("TX Timeout......");
}

2022-12-08

2.1 LoRa Multipunto - HELTEC CubeCell Receptor.ino

Configuración de como receptor de módulo de desarrollo CubeCell

Resultados en monitor serie

RX Paquete "312139" , tamanio 6, Rssi -22 
RX Paquete "312140" , tamanio 6, Rssi -23

Instrucciones

Las instrucciones de dividen en el bloque principal, el procedimiento LoRa Recepcion, separados en cada pestaña del IDE Arduino

Bloque principal

Declara las librerias para el módulo o placa de desarrollo Heltec, se indica los parámetros LoRa como la Banda ISM que para Ecuador es US915, también se establecen las variables para el manejo de los mensajes de recepción, tiempo entre envío de mensajes de número de paquete 'txNumero'

El bucle de configuración setup() inicializa el módulo y el de operación loop() revisa los tiempos en los que se debe realizar el envío del mensaje LoRa.

/* LoRa RECIBE
 * Referencia:  https://github.com/HelTecAutomation/ASR650x-Arduino
 */
#include "LoRaWan_APP.h"
#include "Arduino.h"

#ifndef LoraWan_RGB    // LED placa
#define LoraWan_RGB 0
#endif

// LoRa Parametros
#define RF_FREQUENCY       915E6     // Hz
#define TX_OUTPUT_POWER    14        // dBm
#define LORA_BANDWIDTH     0         // [0: 125 kHz, 1: 250 kHz,
                                     //  2: 500 kHz, 3: Reserved]
#define LORA_SPREADING_FACTOR   7    // [SF7..SF12]
#define LORA_CODINGRATE         1    // [1: 4/5,  2: 4/6,
                                     //  3: 4/7,  4: 4/8]
#define LORA_PREAMBLE_LENGTH    8    // Same for Tx and Rx
#define LORA_SYMBOL_TIMEOUT     0    // Symbols
#define LORA_FIX_LENGTH_PAYLOAD_ON  false
#define LORA_IQ_INVERSION_ON        false
#define RX_TIMEOUT_VALUE            1000
#define BUFFER_SIZE     30      // Define tamanio payload
char txpacket[BUFFER_SIZE];     // cadena de caracteres
char rxpacket[BUFFER_SIZE];
static RadioEvents_t RadioEvents;
int16_t Rssi,rxSize;

void setup() {
  Serial.begin(115200);
  Rssi=0;
  RadioEvents.RxDone = OnRxDone;
  Radio.Init( &RadioEvents );
  Radio.SetChannel( RF_FREQUENCY );
  Radio.SetRxConfig( MODEM_LORA, LORA_BANDWIDTH, LORA_SPREADING_FACTOR,
                    LORA_CODINGRATE, 0, LORA_PREAMBLE_LENGTH,
                    LORA_SYMBOL_TIMEOUT, LORA_FIX_LENGTH_PAYLOAD_ON,
                    0, true, 0, 0, LORA_IQ_INVERSION_ON, true );
  turnOnRGB(COLOR_RECEIVED,0); //change rgb color
  Serial.println("Modo RX Receptor");
}

void loop(){
  Radio.Rx( 0 );
  delay(100);
  Radio.IrqProcess( );
}

Procedimiento de recepción de paquete LoRa

Se desarrollan las instrucciones del receptor en otra pestaña para simplificar el procesamiento del paquete recibido

void OnRxDone( uint8_t *payload, uint16_t size, int16_t rssi, int8_t snr ){
  turnOnRGB(COLOR_RECEIVED,0);
  Rssi = rssi;     // nivel de recepcion
  rxSize = size;
  memcpy(rxpacket, payload, size );
  rxpacket[size] = '\0';  //añade fin de cadena

  Radio.Sleep( );

  Serial.printf("\rRX Paquete \"%s\", tamanio %d, Rssi %d , \r\n",
                rxpacket,rxSize,Rssi);
}

2022-12-08

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

2022-08-19

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

2022-04-17

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

2021-12-23

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

2021-12-03

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)

2021-12-03

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

2021-11-05

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

2021-11-05