Edison Del Rosario

2022-12-08

5. LoRaWan – HELTEC CubeCell Rssi PracticaLab.ino

Ejercicio de conexión a LoRaWan con ChirpStack y HomeAssistant

Bloque principal

// Lectura de Rssi Snr, datarate Up/Downlink
// Datos Downlink de la trama de confirmación anterior
// http://blog.espol.edu.ec/girni/lorawan-enlaces-up-down-archivo-ino/
#include "LoRaWan_APP.h"
#include "Arduino.h"

/* set LoraWan_RGB to Active,the RGB active in loraWan
 * red   |sending;   purple | joined done;
 * blue  |RxWindow1; yellow | means RxWindow2;
 * green | received done;
 */
/* Conexión LoRa: OTAA parametros*/
uint8_t devEui[] = { 0xa5, 0x3e, 0xc6, 0x15, 0xae, 0xde, 0x3f, 0x00 };
uint8_t appEui[] = { 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00 };
uint8_t appKey[] = { 0x88, 0xbe, 0x25, 0xca, 0x2c, 0xcf, 0x31, 0x85,
                     0x51, 0x2d, 0xee, 0xe2, 0x80, 0x31, 0x8e, 0x01 };
/* ABP parametros*/
uint8_t nwkSKey[] = { 0x15, 0xb1, 0xd0, 0xef, 0xa4, 0x63, 0xdf, 0xbe,
                      0x3d, 0x11, 0x18, 0x1e, 0x1e, 0xc7, 0xda,0x85 };
uint8_t appSKey[] = { 0x47, 0xdc, 0xac, 0x5f, 0xc2, 0x32, 0x24, 0x31, 
                      0xdf, 0xf1, 0xff, 0xf9, 0x46, 0xe5, 0x2e, 0x17 };
uint32_t devAddr =  ( uint32_t )0x007bc4200;
/*LoraWan channelsmask, default channels 0-7*/ 
uint16_t userChannelsMask[6]={ 0x00FF,0x0000,0x0000,0x0000,0x0000,0x0000 };

/*Select in arduino IDE tools*/
LoRaMacRegion_t loraWanRegion = ACTIVE_REGION;
DeviceClass_t  loraWanClass = LORAWAN_CLASS;
bool overTheAirActivation = LORAWAN_NETMODE;
bool loraWanAdr = LORAWAN_ADR;
bool keepNet = LORAWAN_NET_RESERVE;
bool isTxConfirmed = LORAWAN_UPLINKMODE;

uint8_t duermemin = 0; //15
uint8_t duermeseg = 30; //0

uint32_t appTxDutyCycle = (duermemin*60 + duermeseg)*1000; // min*seg*ms

uint8_t appPort = 4; /* Application port */
/* trials to transmit frame, if didn't receive ack.
 * The MAC performs a datarate adaptation,
 * Tx nb|Data Rate
 * -----|----------
 * 1    |DR           * 5    | max(DR-2,0)
 * 2    |DR           * 6    | max(DR-2,0)
 * 3    |max(DR-1,0)  * 7    | max(DR-3,0)
 * 4    |max(DR-1,0)  * 8    | max(DR-3,0)
*/
uint8_t confirmedNbTrials = 4;

// Ack parametros de recepción
uint8_t confirmaRssi = 0;
uint8_t confirmaSnr = 0;
uint8_t confirmaDatarate = 0;

uint8_t itera = 0;
uint8_t estado = 0x00; //0x00, 0x01,"OFF","ON"

void setup() {
	Serial.begin(115200);
 
#if(AT_SUPPORT)
	enableAt();
#endif

  // OLED display status
  //LoRaWAN.displayMcuInit();
  
	deviceState = DEVICE_STATE_INIT;
	//LoRaWAN.ifskipjoin(); //if joinned,skip
}

void loop() {
  Serial.print(".");
  itera = itera + 1;
  if (itera>6){
    itera = 0;
    Serial.println(" ");
  }
	switch( deviceState ) {
		case DEVICE_STATE_INIT: {
#if(LORAWAN_DEVEUI_AUTO)
			LoRaWAN.generateDeveuiByChipID();
#endif
#if(AT_SUPPORT)
			getDevParam();
#endif
			printDevParam();
			LoRaWAN.init(loraWanClass,loraWanRegion);
			deviceState = DEVICE_STATE_JOIN;
			break;
		}
		case DEVICE_STATE_JOIN: {
      //LoRaWAN.displayJoining();
			LoRaWAN.join();
			break;
		}
		case DEVICE_STATE_SEND:	{
      //LoRaWAN.displaySending();
			prepareTxFrame( appPort );
			LoRaWAN.send();
			deviceState = DEVICE_STATE_CYCLE;
			break;
		}
		case DEVICE_STATE_CYCLE: {
			// Schedule next packet transmission
			txDutyCycleTime = appTxDutyCycle + randr( 0, APP_TX_DUTYCYCLE_RND );
			LoRaWAN.cycle(txDutyCycleTime);
			deviceState = DEVICE_STATE_SLEEP;
			break;
		}
		case DEVICE_STATE_SLEEP: {
      //LoRaWAN.displayAck();
			LoRaWAN.sleep();
			break;
		}
		default: {
			deviceState = DEVICE_STATE_INIT;
			break;
		}
	}
}

LoraWan transmite

/* Prepares the payload of the frame */
static void prepareTxFrame( uint8_t port ) {
  // enciende sensor
  pinMode(Vext, OUTPUT);
  digitalWrite(Vext, LOW);
  
  //Lectura de Sensor

  // apaga sensor
  digitalWrite(Vext, HIGH);
  
  // lectura de bateria  
  uint16_t batteryVoltage = getBatteryVoltage();
  unsigned char *puc;

  // trama
  appDataSize = 5;
  appData[0] = confirmaRssi; //Ack leido en dispositivo
  appData[1] = confirmaSnr;
  appData[2] = confirmaDatarate;
  appData[3] = (uint8_t)batteryVoltage;
  appData[4] = (uint8_t)(batteryVoltage>>8);

  Serial.print("%, Bateria = ");
  Serial.println(batteryVoltage);
}

LoRWan Recibe

//downlink data handle function example
void downLinkDataHandle(McpsIndication_t *mcpsIndication) {
  // revisa parametros
  // Serial.print("\nLlegó un mensaje para dispositivo...");
  // Serial.print("Rssi: ");
  // Serial.println(mcpsIndication->Rssi);

  // Serial.printf("+REV DATA:%s,RXSIZE %d,PORT %d\r\n",
  //               mcpsIndication->RxSlot ? "RXWIN2" : "RXWIN1",
  //               mcpsIndication->BufferSize, mcpsIndication->Port);
  // Serial.print("+REV DATA:");
  // for (uint8_t i = 0; i < mcpsIndication->BufferSize; i++) {
  //   Serial.printf("%02X", mcpsIndication->Buffer[i]);
  // }

  estado = uint8_t(mcpsIndication->Buffer[0]);
  Serial.print("uplink: rssi = -");
  Serial.println(estado);
  Serial.println();
//  uint32_t color = mcpsIndication->Buffer[0] << 16 | mcpsIndication->Buffer[1] << 8 | mcpsIndication->Buffer[2];
//#if (LoraWan_RGB == 1)
//  turnOnRGB(color, 5000);
//  turnOffRGB();
//#endif
}

LoRaWan recibe confirma (Ack)

void downLinkAckHandle(McpsIndication_t *mcpsIndication){
  confirmaRssi = uint8_t(abs(mcpsIndication->Rssi));
  confirmaSnr  = uint8_t(mcpsIndication->Snr);
  confirmaDatarate = uint8_t(mcpsIndication->RxDoneDatarate);
  //Serial.println(' ');
  //Serial.print(" ack received(rssi,snd,datarate): -");
  //Serial.print(confirmaRssi);Serial.print(" ,");
  //Serial.print(confirmaSnr);Serial.print(" ,");
  //Serial.println(confirmaDatarate);
}

2022-12-082022-12-08

4. LoRa Multipunto – HELTEC CubeCell SemiDuplex.ino con direcciones

Aunque se reciben todos los mensajes en el radio, se requiere discriminar los mensajes que son dirigidos hacia el dispositivo local.

Se modifica las intrucciones de recepción añadiendo un bloque para discriminar si el mensaje es para el dispositivo local. Esto implica descomponer la trama enviada en sus partes y convertir al tipo de dato a usar.

Las operaciones se realizan por caracter o byte con el objetivo de establecer el mecanismo a usar cuando se usa LoRaWan

Resultados obtenidos en puerto serial

TX Paquete "312159" , tamano 6 bytes , proximo en 9341 ms 
    RX Paquete "213169" , tamanio 6, Rssi -20 
    Lee mensaje enviado por: 21 valor mensaje: 69
TX Paquete "312160" , tamano 6 bytes , proximo en 8062 ms 
    RX Paquete "213170" , tamanio 6, Rssi -20 
    Lee mensaje enviado por: 21 valor mensaje: 70

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);
                
  //revisa direcciones
  char envia[4];  byte dir_envia;
  char recibe[4]; byte dir_recibe;
  char msj[4];    byte msj_valor;
  envia[0] = rxpacket[0];
  envia[1] = rxpacket[1];
  envia[2] = '\0';  //añade fin de cadena
  recibe[0] = rxpacket[2];
  recibe[1] = rxpacket[3];
  recibe[2] = '\0';  //añade fin de cadena
  msj[0] = rxpacket[4];
  msj[1] = rxpacket[5];
  msj[2] = '\0';  //añade fin de cadena
  
  //convierte a tipo de datos
  dir_envia  = byte(atoi(envia));
  dir_recibe = byte(atoi(recibe));
  msj_valor  = byte(atoi(msj));
  
  // Muestra o discrimina mensaje en Serial-USB/Pantalla
  if (dir_recibe = dir_local){
    Serial.print("    -- Lee mensaje enviado por: ");
    Serial.print(dir_envia);
    Serial.print(" valor mensaje: ");
    Serial.println(msj_valor);
  }
}

2022-12-082022-12-08

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-082022-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-082022-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-08-192022-08-24

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-04-172022-07-25

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

2021-12-232022-07-26

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-032022-07-25

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-032022-07-25

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)