Proyectos realizados o en desarrollo sobre tecnologías inalámbricas, plataformas tecnológicas.
Categoría: LoRa - Localización Rssi
-
6.1 Puntos de prueba vs Modelo de propagación
Una revisión del modelo de propagación encontrado con los puntos de entrenamiento se puede realizar con puntos de prueba. Las lecturas de los puntos de prueba se realizaron dos semanas después de los puntos de entrenamiento.
Los puntos de prueba tienen una lectura a cada baliza, a diferencia de al menos cien de los puntos de entrenamiento.
En la revisión se encuentra que los puntos de prueba se encuentran en el rango de las mediciones del modelo, sin embargo los puntos se encuentran distribuidos en la gráfica desplazados hacia arriba del modelo, indicando que es necesario incorporar más variables para ajustar el modelo, por ejemplo la temperatura del entorno u otra variable que se determinaría en una siguiente fase del experimento.
Para cuantificar el desplazamiento Δrssi, se calcula el error entre cada punto y la ecuación modelo. Bajo el concepto que los errores se distribuyen equitativamente sobre y debajo la ecuación, se obtiene el valor de la mediana que indica el desplazamiento más optimo como ajuste al modelo bajo las nuevas condiciones.
Revisando la aplicabilidad del modelo, se obtiene el coeficiente de correlación para los puntos de prueba clasificados en gtwFIEC, que es de -0.88 que permite estimar que se puede aplicar relación lineal con pendiente negativa en el intervalo.
gtwFIEC s0 r0 error_mediana: 5.99 correlacion: -0.88 gtwFIEC s0 r1 error_mediana: 7.14 correlacion: -0.68 gtwFIEC s0 r2 error_mediana: 1.98 correlacion: -0.76
Se puede proseguir con los otros puntos de prueba para las otras balizas, encontrando que en la baliza gtwRECT los mejores valores de correlación con los datos disponibles tan dado en el área de vegetación sin sombra. El modelo estimado en vegetación para la zona con sombra debe ser revisados pues hay un valores cercanos a cero en correlación, por lo que la linealización con los datos puede ser insuficiente para estimar que hay una relación lineal entre las dos variables. Se requiere en ese caso expandir la toma de muestra y revisar con mas datos el modelo.
Desplazamientos de puntos de prueba respecto al modelo: [baliza,sector,intervalo,drssi] gtwRECT s0 r0 error_mediana: 7.8 correlacion: -0.54 gtwRECT s0 r1 error_mediana: 5.55 correlacion: -0.13 gtwRECT s0 r2 error_mediana: 9.3 correlacion: -0.47 gtwRECT s1 r0 error_mediana: 3.7 correlacion: 0.13 gtwRECT s1 r1 error_mediana: 3.7 correlacion: 0.13 gtwFIEC s0 r0 error_mediana: 5.99 correlacion: -0.88 gtwFIEC s0 r1 error_mediana: 7.14 correlacion: -0.68 gtwFIEC s0 r2 error_mediana: 1.98 correlacion: -0.76 gtwFCNM s0 r0 error_mediana: 8.56 correlacion: -0.76 gtwFCNM s0 r1 error_mediana: 6.59 correlacion: -0.75 gtwFCNM s0 r2 error_mediana: 9.98 correlacion: -0.29
Algoritmo en Python
Se adjunta el algoritmo correspondiente:
# LoRa-Multipunto, Revisa grafica Rssi vs distancia # con puntos de prueba 'CIRC' # linealización Rssi vs log10(distancia) por mínimos cuadrados # Graficas 2D y 3D # Girni 2020-01-30 propuesta: edelros@espol.edu.ec import numpy as np import pandas as pd import json import matplotlib.pyplot as plt import girni_lora_libreria as girni # INGRESO # archivos de entrada modo = 'rx' medida = 'rssi' descriptor = 'mean' arch_ecuaciones = 'rsmP07_ecuacionSector01.json' arch_medidaubica = 'rsmP06_'+medida+'UbicaCirc01.txt' # archivos de salida arch_ecuaciones2 = 'rsmP07_ecuacionSector02.json' baliza = {'d1':'gtwRECT', 'd2':'gtwFIEC', 'd3':'gtwFCNM'} # Parámetros de grafica tipograf = '2D' # '2D','3D' escala = 'log' # 'normal','log' escalabase = 10 # 10, np.exp() casicero = 1e-4 precision = 2 intersectar = 0 # 0:Falso, 1: Verdadero # Referencias de gráfica grupo = ['FIEC' ,'FCNM' ,'RECT','CIRC'] colores = ['green','orange','grey','magenta'] tipo = ['punto','1m' ,'gtw','dispositivo'] marcas = [ 'o','D' ,'D' ,'*' ] mostrargrpeti = ['FIEC','FCNM','RECT'] mostrartipeti = ['1m','gtw'] # PROCEDIMIENTO # leer datos with open(arch_ecuaciones) as json_file: ecuacion = json.load(json_file) tabla = pd.read_csv(arch_medidaubica, index_col='etiqueta') tabla = pd.DataFrame(tabla) baliza_key = list(baliza.keys()) baliza_val = list(baliza.values()) eq_graf = {} # Revisa un sector_estimado for cualbaliza in baliza: sector_baliza = 'sector_'+cualbaliza tabla[sector_baliza] = 0 for cadapunto in tabla.index: for cualbaliza in baliza: # evalua en sector sector_baliza = 'sector_'+cualbaliza donde = baliza_key.index(cualbaliza) # coordenadas baliza para identificar ángulo b_este = tabla['c_este'][baliza[cualbaliza]] b_norte = tabla['c_norte'][baliza[cualbaliza]] # coordenadas del punto p_este = tabla['c_este'][cadapunto] p_norte = tabla['c_norte'][cadapunto] dx = p_este-b_este dy = p_norte-b_norte theta = np.arctan2(dy,dx) if theta<0: theta = theta + 2*np.pi sectores = ecuacion[baliza[cualbaliza]]['sector_rad'] nsectores = len(sectores) otrosector = '' for i in range(0,nsectores-1,1): a = sectores[i] b = sectores[i+1] if theta>=a and theta<b: otrosector = i+1 if otrosector !='': tabla.loc[cadapunto,'sector_'+cualbaliza] = otrosector # analiza datos hacia una baliza for unabaliza in ecuacion: donde = baliza_val.index(unabaliza) cualbaliza = baliza_key[donde] eq_graf[unabaliza] = {} for unsector in list(ecuacion[unabaliza].keys()): if unsector != 'sector_rad': eq_graf[unabaliza][unsector] = {} for i_eq in list(ecuacion[unabaliza][unsector].keys()): eq_graf[unabaliza][unsector][i_eq] = {} # puntos de tabla for unsector in list(ecuacion[unabaliza].keys()): if unsector != 'sector_rad': xi_p = [] ; yi_p = [] xi = [] ; yi = [] ; etiqueta = [] columna = medida+'_'+modo+'_'+cualbaliza for cadapunto in tabla.index: p_rssi = np.round(tabla[columna][cadapunto], precision) dist = np.round(tabla['dist_'+cualbaliza][cadapunto], precision) ensector = 's'+str(tabla['sector_'+cualbaliza][cadapunto]) LOS = tabla['LOS_'+cualbaliza][cadapunto] if not(np.isnan(p_rssi)) and dist>1 and LOS==1 and unsector==ensector: yi_p.append(p_rssi) xi_p.append(dist) etiqueta.append(cadapunto) if not(dist in xi): xi.append(dist) yi.append(p_rssi) xi_p = np.array(xi_p) yi_p = np.array(yi_p) ordenar = np.argsort(xi_p) xi_p = list(xi_p[ordenar]) yi_p = list(yi_p[ordenar]) etiqueta = np.array(etiqueta) etiqueta = etiqueta[ordenar] xi = np.array(xi) yi = np.array(yi) ordenar = np.argsort(xi) xi = xi[ordenar] yi = yi[ordenar] i_eq = 'r0' # todos los puntos eq_graf[unabaliza][unsector][i_eq]['xi'] = xi_p.copy() eq_graf[unabaliza][unsector][i_eq]['yi'] = yi_p.copy() eq_graf[unabaliza][unsector][i_eq]['etiqueta'] = etiqueta.copy() # ecuacion linealizada for i_eq in list(ecuacion[unabaliza][unsector].keys()): a = ecuacion[unabaliza][unsector][i_eq]['intervalox'][0] b = ecuacion[unabaliza][unsector][i_eq]['intervalox'][1] alpha = ecuacion[unabaliza][unsector][i_eq]['alpha'] beta = ecuacion[unabaliza][unsector][i_eq]['beta'] std = ecuacion[unabaliza][unsector][i_eq]['error_std'] fdist = lambda d: -10*alpha*(np.log10(d))+beta subintervalo = (xi >= a) & (xi <= b) xi_s = xi[subintervalo] yi_s = yi[subintervalo] # correlacion de los puntos correlacion = np.corrcoef(xi_s,yi_s)[0,1] # dibujar la línea en intervalo xi_sub = xi[subintervalo] if not(a in xi_sub): xi_sub = np.concatenate(([a],xi_sub),axis=0) if not(b in xi_sub): xi_sub = np.concatenate((xi_sub,[b]),axis=0) yi_sub = fdist(xi_sub) eq_graf[unabaliza][unsector][i_eq]['xi_graf'] = xi_sub.copy() eq_graf[unabaliza][unsector][i_eq]['yi_graf'] = yi_sub.copy() # Errores en subintervalo mediana = 0 if intersectar == 0: yie_s = yi_s - fdist(xi_s) mediana = np.median(yie_s) else: if i_eq=='r0': yie_s = yi_s - fdist(xi_s) mediana = np.median(yie_s) else: mediana = ecuacion[unabaliza][unsector]['r0']['desplaza'] eq_graf[unabaliza][unsector][i_eq]['desplaza'] = mediana eq_graf[unabaliza][unsector][i_eq]['correlacion'] = correlacion ecuacion[unabaliza][unsector][i_eq]['desplaza'] = mediana # SALIDA print('Desplazamientos de puntos de prueba') print(' respecto al modelo:') print('[baliza,sector,intervalo,drssi]') for unabaliza in ecuacion: for unsector in list(ecuacion[unabaliza].keys()): if unsector != 'sector_rad': for i_eq in list(ecuacion[unabaliza][unsector].keys()): mediana = ecuacion[unabaliza][unsector][i_eq]['desplaza'] mediana = np.round(mediana,precision) correlacion = eq_graf[unabaliza][unsector][i_eq]['correlacion'] correlacion = np.round(correlacion,precision) texto = unabaliza + ' '+unsector + ' ' + i_eq + ' ' texto = texto + ' error_mediana: '+ str(mediana) texto = texto + ' correlacion: '+str(correlacion) print(texto) # salida hacia archivo with open(arch_ecuaciones2, 'w') as outfile: json.dump(ecuacion, outfile) # GRAFICA # Referencias para gráfica grupo = ['FIEC' ,'FCNM' ,'RECT','CIRC'] colores = ['green','orange','grey','magenta'] tipo = ['punto','1m' ,'gtw','dispositivo'] marcas = [ 'o','D' ,'D' ,'*' ] mostrargrpeti = ['FIEC','FCNM','RECT'] mostrartipeti = ['1m','gtw'] if tipograf=='2D': for unabaliza in eq_graf: for unsector in eq_graf[unabaliza]: figura,grafica = plt.subplots() if escala == 'log': grafica.set_xscale(escala,base=escalabase) eq_graf2 = {} # todos los puntos unintervalo = 'r0' xi_p = eq_graf[unabaliza][unsector][unintervalo]['xi'] yi_p = eq_graf[unabaliza][unsector][unintervalo]['yi'] etiqueta = eq_graf[unabaliza][unsector][unintervalo]['etiqueta'] grafica.scatter(xi_p,yi_p,marker='.') m = len(xi_p) for i in range(0,m,1): grafica.annotate(etiqueta[i], (xi_p[i],yi_p[i])) # linea con todos los puntos for i_eq in list(eq_graf[unabaliza][unsector].keys()): fdtxt = ecuacion[unabaliza][unsector][i_eq]['eq_latex'] xi_graf = eq_graf[unabaliza][unsector][i_eq]['xi_graf'] yi_graf = eq_graf[unabaliza][unsector][i_eq]['yi_graf'] a = np.round(ecuacion[unabaliza][unsector][i_eq]['intervalox'][0],precision) b = np.round(ecuacion[unabaliza][unsector][i_eq]['intervalox'][1],precision) desplaza = eq_graf[unabaliza][unsector][i_eq]['desplaza'] desplaza = np.round(desplaza,precision) estilo = 'solid' if i_eq =='r0': estilo = 'dashed' eq_texto = i_eq+': '+fdtxt+' ; ['+str(a)+','+str(b)+']' grafica.plot(xi_graf,yi_graf, label = eq_texto, linestyle = estilo) # para linea desplazada eq_graf2[i_eq] = {'xi':xi_graf, 'yi':yi_graf + desplaza, 'desplaza': desplaza} # lineas de frontera grafica.axvline(a, color='lightblue') valor_frontera = str(np.round(a,precision)) grafica.annotate(valor_frontera, (a,np.max(yi_p)), color='lightblue') grafica.axvline(b, color='lightblue') valor_frontera = str(np.round(b,precision)) grafica.annotate(valor_frontera, (b,np.max(yi_p)), color='lightblue') for i_eq in list(ecuacion[unabaliza][unsector].keys()): estilo = 'solid' if i_eq =='r0': estilo = 'dashed' xi = eq_graf2[i_eq]['xi'] yi = eq_graf2[i_eq]['yi'] desplaza = eq_graf2[i_eq]['desplaza'] texto = i_eq + ' + $\Delta $Rssi' texto = texto + ' ; $ \Delta rssi$= '+ str(desplaza) grafica.plot(xi,yi, label = texto, linestyle = estilo) # etiquetas y títulos grafica.legend() grafica.set_ylabel(medida+'_'+modo) grafica.set_xlabel('distancia') grafica.grid(True,linestyle='dotted', axis='x', which='both') untitulo = unabaliza+' '+unsector+': Puntos de prueba desplazados' grafica.set_title(untitulo) plt.show()
-
3.6 Procesa Datos. Un punto - revisa descriptores estadísticos
El comportamiento de la señal de cada baliza hacia un punto en particular se observa al calcular algunos valores de estadística.
El valor principal para el modelo de localización usado es el promedio de rssi de cada baliza en el punto.
Otro valor comprementario la desviación estándar (std, σ) que es una medida de dispersión alrededor de la media. Estas variaciones estan asociadas a múltiples factores que depende de las características del entorno.
Si hay valores de promedios de un punto que son atípicos a la tendencia entre Rssi y la distancia se puede recurrir a ésta sección para observar el punto.
Con el algoritmo se obtienen los siguientes valores para el punto del ejemplo:
baliza: d1 gtwRECT muestras: 137 tramo: 17 count mean std min 25% 50% 75% max 0 118.0 -128.211864 1.960515 -132.0 -130.0 -128.5 -127.0 -123.0 baliza: d2 gtwFIEC muestras: 144 tramo: 15 count mean std min 25% 50% 75% max 0 118.0 -128.211864 1.960515 -132.0 -130.0 -128.5 -127.0 -123.0 baliza: d3 gtwFCNM muestras: 118 tramo: 9 count mean std min 25% 50% 75% max 0 118.0 -128.211864 1.960515 -132.0 -130.0 -128.5 -127.0 -123.0
Datos de ingreso
El primer valor que se requiere es el nombre del punto, ejemplo 'FIEC102', luego el modo de lectura (rx,tx') y la medida (rssi, snr).
Los datos se toman del archivo que integra en una tbala los datos del punto, en formato json del proceso anterior.
Se usa el diccionario que identifica a las balizas para realizar la tabulación de los descriptores.
# INGRESO cadapunto = 'FIEC102' modo = 'rx' medida = 'rssi' descriptor = 'mean' # Archivos de datos arch_detalles = 'rsmP02_tabla_puntosdatos01.json' # Referencias baliza = {'d1':'gtwRECT', 'd2':'gtwFIEC', 'd3':'gtwFCNM'} tolera = 1e-8
Los resultados se complementan con la gráfica que muestra la distribución de probabilidad para cada baliza, marcando también el valor de promedio y la suma de promedio mas desvición estándar.
El análisis del comportamiento de la señal de cada punto tiene relación con el entorno del punto dominado por edificios o vegetación.
Algoritmo en Python
# Revisar un punto, descriptores estadísticos # requiere un punto, modo, medida # nombre del archivo de la tabla detalle # Girni 2020-10-07 edelros@espol.edu.ec import numpy as np import pandas as pd import json import matplotlib.pyplot as plt # INGRESO cadapunto = 'FIEC102' modo = 'rx' medida = 'rssi' descriptor = 'mean' # Archivos de datos arch_detalles = 'rsmP02_tabla_puntosdatos01.json' # Referencias baliza = {'d1':'gtwRECT', 'd2':'gtwFIEC', 'd3':'gtwFCNM'} tolera = 1e-8 # PROCEDIMIENTO baliza_key = list(baliza.keys()) baliza_val = list(baliza.values()) # leer datos with open(arch_detalles) as json_file: puntodato = json.load(json_file) punto = {} for cadabaliza in baliza: punto[cadabaliza] = {} # revisa un punto unpunto = puntodato[modo][cadapunto][cadabaliza][medida+'_'+modo] unpunto = np.array(unpunto) ordenar = np.argsort(unpunto) unpunto = unpunto[ordenar] puntopd = pd.DataFrame(unpunto) describe = puntopd.describe() describe = describe.T pmin = describe['min'][0] pmax = describe['max'][0] tramo = int(pmax-pmin) conteo = np.zeros(tramo+1,dtype=int) posicion = np.arange(pmin,pmax+1,1) posicion = list(posicion) conteo = list(conteo) for valor in unpunto: donde = posicion.index(valor) conteo[donde] = conteo[donde] + 1 punto[cadabaliza] = {'muestras': len(unpunto), 'tramo': tramo, 'posicion': posicion.copy(), 'conteo': conteo.copy(), 'describe': describe} # SALIDA for cadabaliza in baliza: print('baliza:', cadabaliza, baliza[cadabaliza]) print('muestras: ', punto[cadabaliza]['muestras'], 'tramo: ' , punto[cadabaliza]['tramo']) print(describe) # grafica precision = 2 figura,grafica = plt.subplots() n_baliza = len(baliza) for cadabaliza in baliza: p_media = np.round(punto[cadabaliza]['describe']['mean'][0],precision) p_std = np.round(punto[cadabaliza]['describe']['std'][0],precision) muestras = punto[cadabaliza]['muestras'] posicion = punto[cadabaliza]['posicion'] conteo = np.array(punto[cadabaliza]['conteo'])/muestras conteomax = np.max(conteo) #grafica.plot(posicion,conteo,'o',label = cadabaliza) grafica.plot(posicion,conteo,'-',label = cadabaliza) grafica.axvline(p_media, color = 'green', linestyle = 'dashed') texto = cadabaliza + ': ' + str(p_media) grafica.annotate(cadabaliza+': '+str(p_media),(p_media,conteomax), color='green') grafica.axvline(p_media+ p_std, color ='orange', linestyle = 'dotted') texto = texto + '+' + str(p_std) grafica.annotate(texto,(p_media+p_std,conteomax-0.015), color='orange') grafica.axvline(p_media- p_std, color ='orange', linestyle = 'dotted') grafica.axvline(p_media+ 2*p_std, color ='yellow', linestyle = 'dotted') grafica.axvline(p_media- 2*p_std, color ='yellow', linestyle = 'dotted') grafica.set_xlabel('rssi') grafica.set_ylabel('pmf') texto = 'Distribución de probabilidad: '+cadapunto+' '+medida grafica.set_title(texto) grafica.grid(True,linestyle='dotted', axis='y', which='both') grafica.legend() plt.show()
-
5.5 Localiza por Trilateración – Revisa un punto
Para revisar los resultados de un solo punto, de todos los resultados anteriores, se realiza una gráfica que muestre los círculos de cada baliza con los radios de las distancias estimadas.
La gráfica permite observar los detalles de trilateración para ese punto en particular.
Resultados: baricentro: [ 614804.97 9762762.46] bar_error: 72.84 poligono x: [614760.66 614776.52 614877.73] poligono y: [9762742.69 9762778.91 9762765.77]
Por ejemplo para el caso del punto 'FIEC130' del que no se registró señal desde 'gtwFCNM' se tiene:
Algoritmo Python
# trilaterando - gráfica de un punto # Girni 2020-10-07 propuesta: edelros@espol.edu.ec import numpy as np import matplotlib.pyplot as plt import pandas as pd import girni_lora_libreria as girni # INGRESO cadapunto = 'FIEC111' modo = 'rx' medida = 'rssi' arch_trilatera = 'rsmP07_'+medida+'trilateraSector01.txt' # Referencias baliza = {'d1':'gtwRECT', 'd2':'gtwFIEC', 'd3':'gtwFCNM'} tolera = 1e-8 # PROCEDIMIENTO baliza_key = list(baliza.keys()) baliza_val = list(baliza.values()) # leer datos tabla = pd.read_csv(arch_trilatera, index_col='etiqueta') tabla = pd.DataFrame(tabla) # Revisa indices existe = 0 if cadapunto in tabla.index: # coordenadas de baliza coordbaliza = {} for cadabaliza in baliza: cualbaliza = baliza[cadabaliza] coordbaliza[cadabaliza] = [tabla['c_este'][cualbaliza]] coordbaliza[cadabaliza].append(tabla['c_norte'][cualbaliza]) existe = 1 # coordenadas GPS del punto c_este = tabla['c_este'][cadapunto] c_norte = tabla['c_norte'][cadapunto] # coordenadas por trilateración punto = {} ubicados = np.zeros(4,dtype=int) for i in range(0,3,1): punto[i] = {} cuenta = 0 for unabaliza in baliza: punto[i][unabaliza] = tabla['dist'+str(i)+'_'+unabaliza][cadapunto] # si existe distancia if not(np.isnan(punto[i][unabaliza])): cuenta = cuenta + 1 if cuenta>0: ubicados[3]= ubicados[3] + 1 encontrado = tabla['encontrado'][cadapunto] if not(np.isnan(encontrado)): encontrado = int(encontrado) else: encontrado = 0 # trilatera el punto raiztodas = girni.intersectacirculos(punto[encontrado],coordbaliza,tolera) todasx = raiztodas[0] todasy = raiztodas[1] if len(todasx)>0: resultado = girni.trilatera(punto[encontrado],coordbaliza,tolera) # SALIDA precision=2 np.set_printoptions(precision) if existe ==1: if len(todasx)>0: print('intersectados: ') print(np.array(todasx)) print(np.array(todasy)) print('\nResultados: ') print('baricentro: ',np.array(resultado['baricentro'])) print('bar_error: ',np.round(resultado['barerror'],precision)) print('poligono x: ',np.array(resultado['poligono'][0])) print('poligono y: ',np.array(resultado['poligono'][1])) else: print('El punto indicado no existe') # Grafica if existe ==1: figura = plt.figure() grafica = figura.add_subplot(111) grafica.set_aspect('equal', adjustable='box') # punto con GPS texto = '['+str(np.round(c_este,2))+', '+str(np.round(c_norte,2))+']' grafica.scatter(c_este,c_norte, marker='D', color = 'grey', label = texto) grafica.annotate('gps',(c_este,c_norte), color ='grey') # cotas de gráfica unalista = list(coordbaliza.values()) unalista = np.array(unalista) cotax = [np.min(unalista[:,0]), np.max(unalista[:,0])] cotay = [np.min(unalista[:,1]), np.max(unalista[:,1])] cotadx = 0.15*(cotax[1]-cotax[0]) cotady = 0.15*(cotay[1]-cotay[0]) grafica.set_xlim(cotax[0]-cotadx,cotax[1]+cotadx) grafica.set_ylim(cotay[0]-cotady,cotay[1]+cotady) # raices k = len(todasx) if k>0: grafica.scatter(todasx,todasy, marker='+', color = 'magenta') poligono = resultado['poligono'] p = len(poligono[0]) for j in range(0,p,1): grafica.scatter(poligono[0][j],poligono[1][j], marker='*', color = 'blue') # coodenada con trilatera baricentro = resultado['baricentro'] sumabaricentro = np.sum(baricentro) if not(np.isnan(sumabaricentro)): texto = '['+str(np.round(baricentro[0],2)) texto = texto +', '+str(np.round(baricentro[1],2))+']' grafica.scatter(baricentro[0],baricentro[1], marker='D', color = 'red', label = texto) grafica.annotate('trilatera', (baricentro[0],baricentro[1]), color ='red') uncirculo = plt.Circle((baricentro[0],baricentro[1]), resultado['barerror'], ec = 'red', fc='None', linestyle ='dashed') grafica.add_artist(uncirculo) # circulos colores = ['blue','orange','lightgreen','grey'] j = 0 for fila in coordbaliza: xc = coordbaliza[fila][0] yc = coordbaliza[fila][1] r = punto[encontrado][fila] uncolor = colores[j] grafica.scatter(xc,yc, marker = '+', color = uncolor) grafica.annotate(fila,(xc,yc)) uncirculo = plt.Circle((xc,yc),r, ec = uncolor, fc='None', linestyle ='dotted') grafica.add_artist(uncirculo) j = j + 1 plt.title('Trilatera: ' + cadapunto) plt.legend() plt.xlabel('coordenada este') plt.ylabel('coordenada norte') plt.show()
-
5.4 Localiza por Trilateración – Sectores e Intervalos. Algoritmo Python
Los sectores en una baliza se incorporan en el algoritmo como la estimación de ubicación del punto usando las otras balizas.
Para un punto de ejemplo, dado que baliza 'gtwRECT' tiene sectores, se usan las otras balizas 'gtwFIEC' y 'gtwFCNM' para estimar la ubicación y el ángulo que permite aproximar el sector donde se ubicaría el punto.
Este sector se usa para aplicar la ecuación que corresponde al 'gtwFIEC'.
El desarrollo continua desde la sección anterior con solo intervalos y se incorpora el concepto.
Algoritmo Python
# Localización por trilateración # Ubicación de puntos en el mapa # Girni 2020-10-07 propuesta: edelros@espol.edu.ec import numpy as np import sympy as sym import matplotlib.pyplot as plt import pandas as pd import json import girni_lora_libreria as girni # INGRESO # revisar parametros al inicio modo = 'rx' medida = 'rssi' # archivos de entrada arch_medUbAtrib = 'rsmP06_'+medida+'UbicaUsarSector01.txt' arch_ecuaciones = 'rsmP06_ecuacionSector01.json' # archivos de salida arch_trilatera = 'rsmP07_'+medida+'trilateraSector01.txt' # Parametros de grafica mostrargrp = ['FIEC'] #['FIEC','FCNM','RECT','CIRC'] mostrartip = ['punto'] mostrarfuera = [0,1] #[0,1] tolera_error = 2 # decimales en error # Referencias baliza = {'d1':'gtwRECT', 'd2':'gtwFIEC', 'd3':'gtwFCNM'} # Parámetros de grafica tipograf = '2D' # '2D','' sin grafica # PROCEDIMIENTO baliza_key = list(baliza.keys()) baliza_val = list(baliza.values()) # leer datos with open(arch_ecuaciones) as json_file: ecuacion = json.load(json_file) tabla = pd.read_csv(arch_medUbAtrib, index_col='etiqueta') tabla = pd.DataFrame(tabla) # coordenadas de baliza coordbaliza = {} for cadabaliza in baliza: cualbaliza = baliza[cadabaliza] coordbaliza[cadabaliza] = [tabla['c_este'][cualbaliza]] coordbaliza[cadabaliza].append(tabla['c_norte'][cualbaliza]) # inicializa ubicados,encontrados, intervalos fuera, errores # ubicados tiene [distancia, distancia+1s, distancia +2s,total] ubicados = np.zeros(shape=(3,4),dtype=int) tabla['interv_fuera'] = '' tabla['encontrado'] = np.nan tabla['ensector'] = '' # suma de errores sumatrilat = np.zeros(shape=(2,4),dtype=float) sumagps = np.zeros(shape=(2,4),dtype=float) sumagpsdx = np.zeros(shape=(2,4),dtype=float) sumagpsdy = np.zeros(shape=(2,4),dtype=float) # distancias y errores a cada baliza for unabaliza in baliza_val: donde = baliza_val.index(unabaliza) cualbaliza = baliza_key[donde] columna = medida+'_'+modo+'_'+cualbaliza tabla['dist0_'+cualbaliza] = np.nan tabla['dist1_'+cualbaliza] = np.nan tabla['dist2_'+cualbaliza] = np.nan # distancia estimada por punto for cadapunto in tabla.index: p_rssi = tabla[columna][cadapunto] esgrupo = tabla['grupo'][cadapunto] estipo = tabla['tipo'][cadapunto] unsector = 's0' if not(np.isnan(p_rssi)): ecuacion_rssi = ecuacion[unabaliza][unsector] if (esgrupo in mostrargrp) and (estipo in mostrartip): ubicados[2,donde] = ubicados[2,donde] + 1 dist = girni.dist_rssi(p_rssi,ecuacion_rssi) [distancia,e_mean,e_1std,e_2std,interv_fuera] = dist dist1 = distancia + e_1std dist2 = distancia + e_2std tabla.loc[cadapunto,'dist0_'+cualbaliza] = distancia tabla.loc[cadapunto,'dist1_'+cualbaliza] = dist1 tabla.loc[cadapunto,'dist2_'+cualbaliza] = dist2 if interv_fuera>0: if len(tabla['interv_fuera'][cadapunto])==0: separador='' else: separador=',' tabla.loc[cadapunto,'interv_fuera'] = tabla['interv_fuera'][cadapunto]+separador+cualbaliza procesados = np.max(ubicados[2,:]) # Revisa un sector_estimado columnas = list(tabla.keys()) for cadapunto in tabla.index: for cualbaliza in baliza: sector_baliza = 'sector_'+cualbaliza unsector = 0 if sector_baliza in columnas: unsector = tabla['sector_'+cualbaliza][cadapunto] if unsector!=0: # evalua en sector donde = baliza_key.index(cualbaliza) otrasbalizas = baliza_key.copy() otrasbalizas.pop(donde) # coordenadas de baliza coord_otrasbalizas = {} for cadaotrabaliza in otrasbalizas: cualotrabaliza = baliza[cadaotrabaliza] coord_otrasbalizas[cadaotrabaliza] = [tabla['c_este'][cualotrabaliza]] coord_otrasbalizas[cadaotrabaliza].append(tabla['c_norte'][cualotrabaliza]) # cada punto a i*std punto = {} sumabaricentro = np.nan for i in range(0,3,1): punto[i] = {} cuenta = 0 for cadaotrabaliza in otrasbalizas: punto[i][cadaotrabaliza] = tabla['dist'+str(i)+'_'+cadaotrabaliza][cadapunto] # si existe distancia if not(np.isnan(punto[i][cadaotrabaliza])): cuenta = cuenta + 1 # Hay 3 mediciones a baliza if cuenta==2: localiza = girni.trilatera(punto[i], coord_otrasbalizas, tolera = 10e-4) baricentro = localiza['baricentro'] barerror = localiza['barerror'] poligono = localiza['poligono'] sumabaricentro = np.sum(baricentro) if not(np.isnan(sumabaricentro)): break if not(np.isnan(sumabaricentro)): # coordenadas baliza para identificar ángulo b_este = tabla['c_este'][baliza[cualbaliza]] b_norte = tabla['c_norte'][baliza[cualbaliza]] # coordenadas del punto entre otras balizas p_este = baricentro[0] p_norte = baricentro[1] dx = p_este-b_este dy = p_norte-b_norte theta = np.arctan2(dy,dx) if theta<0: theta = theta + 2*np.pi sectores = ecuacion[baliza[cualbaliza]]['sector_rad'] nsectores = len(sectores) otrosector = '' for i in range(0,nsectores-1,1): a = sectores[i] b = sectores[i+1] if theta>=a and theta<b: otrosector = i+1 if otrosector !='': tabla.loc[cadapunto,'sector_'+cualbaliza] = otrosector # recalcula distancia con sector columna = medida+'_'+modo+'_'+cualbaliza p_rssi = tabla[columna][cadapunto] ecuacion_rssi = ecuacion[baliza[cualbaliza]]['s'+str(otrosector)] dist = girni.dist_rssi(p_rssi,ecuacion_rssi) [distancia,e_mean,e_1std,e_2std,interv_fuera] = dist dist1 = distancia + e_1std dist2 = distancia + e_2std tabla.loc[cadapunto,'dist0_'+cualbaliza] = distancia tabla.loc[cadapunto,'dist1_'+cualbaliza] = dist1 tabla.loc[cadapunto,'dist2_'+cualbaliza] = dist2 # trilateración. localiza cada punto # cada punto a i*std for i in range(0,3,1): tabla['trilat_este_'+str(i)] = np.nan tabla['trilat_norte_'+str(i)] = np.nan tabla['trilat_error_'+str(i)] = np.nan tabla['ubicado_'+str(i)] = np.nan tabla['trilat_gps_error_'+str(i)] = np.nan tabla['trilat_gps_dx_'+str(i)] = np.nan tabla['trilat_gps_dy_'+str(i)] = np.nan # procesa cada punto for cadapunto in tabla.index: esgrupo = tabla['grupo'][cadapunto] # coordenadas GPS del punto c_este = tabla['c_este'][cadapunto] c_norte = tabla['c_norte'][cadapunto] if esgrupo in mostrargrp: # cada punto a i*std punto = {} for i in range(0,3,1): punto[i] = {} cuenta = 0 for unabaliza in baliza: punto[i][unabaliza] = tabla['dist'+str(i)+'_'+unabaliza][cadapunto] # si existe distancia if not(np.isnan(punto[i][unabaliza])): cuenta = cuenta + 1 # Hay 3 mediciones a baliza if cuenta==3: localiza = girni.trilatera(punto[i], coordbaliza, tolera = 10e-4) baricentro = localiza['baricentro'] barerror = localiza['barerror'] poligono = localiza['poligono'] sumabaricentro = np.sum(baricentro) # Hay coordenadas de baricentro if not(np.isnan(sumabaricentro)): tabla.loc[cadapunto,'trilat_este_'+str(i)] = baricentro[0] tabla.loc[cadapunto,'trilat_norte_'+str(i)] = baricentro[1] tabla.loc[cadapunto,'trilat_error_'+str(i)] = np.round(barerror,tolera_error) # error trilatera hacia gps dx = baricentro[0] - c_este dy = baricentro[1] - c_norte error_gps = np.sqrt(dx**2+dy**2) tabla.loc[cadapunto,'trilat_gps_error_'+str(i)] = np.round(error_gps,tolera_error) tabla.loc[cadapunto,'trilat_gps_dx_'+str(i)] = np.round(dx,tolera_error) tabla.loc[cadapunto,'trilat_gps_dy_'+str(i)] = np.round(dy,tolera_error) # cuenta error si encuentra primera vez cond1 = np.isnan(tabla['encontrado'][cadapunto]) cond2 = len(tabla['interv_fuera'][cadapunto])>0 k = int(cond2) #intervalo extendido if cond1: tabla.loc[cadapunto,'encontrado'] = i ubicados[k,i] = ubicados[k,i]+1 sumatrilat[k,i] = sumatrilat[k,i] + np.round(barerror,tolera_error) sumagps[k,i] = sumagps[k,i] + np.round(error_gps,tolera_error) sumagpsdx[k,i] = sumagpsdx[k,i] + np.round(np.abs(dx),tolera_error) sumagpsdy[k,i] = sumagpsdy[k,i] + np.round(np.abs(dy),tolera_error) # contabiliza errores de localizados ubicasuma = np.zeros(2) ubicasumporc = np.zeros(2) ubicaporc = np.zeros(shape=(2,3)) for k in range(0,2,1): ubicados[k,3] = int(np.sum(ubicados[k,0:3])) ubicasuma[k] = np.sum(ubicados[k,0:3]) ubicasumporc[k] = np.round(100*ubicasuma[k]/ubicados[k,3],1) for i in range(0,3,1): ubicaporc[k,i] = np.round(100*ubicados[k,i]/ubicasuma[k],1) if ubicados[k,i]>0: sumatrilat[k,i] = np.round(sumatrilat[k,i]/ubicados[k,i],1) sumagps[k,i] = np.round(sumagps[k,i]/ubicados[k,i],1) sumagpsdx[k,i] = np.round(sumagpsdx[k,i]/ubicados[k,i],1) sumagpsdy[k,i] = np.round(sumagpsdy[k,i]/ubicados[k,i],1) # SALIDA print('Errores localizacion') import prettytable as ptt print('Errores estimado: Cota-Trilatera-polígono y Trilatera_vs_GPS') mostrar = ptt.PrettyTable(['punto','i*std','fuera','trilat', 'gps','gps_dx','gps_dy','sect_d1','sect_d2','sect_d3']) for cadapunto in tabla.index: # selecciona puntos a mostrar encontrado = tabla['encontrado'][cadapunto] esgrupo = tabla['grupo'][cadapunto] interv_fuera = int(len(tabla['interv_fuera'][cadapunto])!=0) cond1 = not(np.isnan(encontrado)) cond2 = esgrupo in mostrargrp cond3 = interv_fuera in mostrarfuera if cond1 and cond2 and cond3: encontrado = int(encontrado) mostrar.add_row([cadapunto, str(int(tabla['encontrado'][cadapunto])), tabla['interv_fuera'][cadapunto], str(tabla['trilat_error_'+str(encontrado)][cadapunto]), str(tabla['trilat_gps_error_'+str(encontrado)][cadapunto]), str(tabla['trilat_gps_dx_'+str(encontrado)][cadapunto]), str(tabla['trilat_gps_dy_'+str(encontrado)][cadapunto]), tabla['sector_d1'][cadapunto], tabla['sector_d2'][cadapunto], tabla['sector_d3'][cadapunto] ]) print(mostrar) print('puntos con medidas: ',procesados) loc_alg = np.round(100*np.sum(ubicados[:,3])/procesados,1) print('localizados con algoritmo: ', np.sum(ubicados[:,3]), ' , '+str(loc_alg)+'%' ) print('Errores Promedio:') resumen = ptt.PrettyTable(['error+i*std','cant','%','trilat', 'gps','gps_dx','gps_dy']) for k in range(0,2,1): if k == 0: texto = 'dentro' if k == 1: texto = 'fuera' resumen.add_row(['interv_'+texto,ubicados[k,3],'', '','','','']) for i in range(0,3,1): resumen.add_row(['error_'+str(i),ubicados[k,i], str(ubicaporc[k,i])+'%', sumatrilat[k,i],sumagps[k,i], sumagpsdx[k,i],sumagpsdy[k,i], ]) print(resumen) # salida hacia archivo tabla.to_csv(arch_trilatera) # Grafica ubicados ------------ # Referencias para gráfica grupo = ['FIEC' ,'FCNM' ,'RECT','CIRC'] colores = ['green','orange','grey','magenta'] tipo = ['punto','1m' ,'gtw','dispositivo'] marcas = [ 'o','D' ,'D' ,'*' ] colorstd = ['lightblue', 'lightgreen','orange'] colorlin = ['lightblue', 'lightgreen','orange'] if tipograf == '2D': figura,grafica = plt.subplots() # balizas for unabaliza in coordbaliza: g_este = coordbaliza[unabaliza][0] g_norte = coordbaliza[unabaliza][1] grafica.scatter(g_este,g_norte, color = 'red', marker = 'D', label = cadapunto) grafica.annotate(unabaliza, (g_este,g_norte)) # Puntos for cadapunto in tabla.index: g_este = tabla['c_este'][cadapunto] g_norte = tabla['c_norte'][cadapunto] # selecciona puntos a mostrar encontrado = tabla['encontrado'][cadapunto] esgrupo = tabla['grupo'][cadapunto] interv_fuera = int(len(tabla['interv_fuera'][cadapunto])!=0) cond1 = not(np.isnan(encontrado)) cond2 = esgrupo in mostrargrp cond3 = interv_fuera in mostrarfuera if cond1 and cond2 and cond3: encontrado = int(encontrado) p_este = tabla['trilat_este_'+str(encontrado)][cadapunto] p_norte = tabla['trilat_norte_'+str(encontrado)][cadapunto] grafica.scatter(p_este,p_norte, color = colorstd[encontrado], label = cadapunto) grafica.plot([p_este,g_este], [p_norte,g_norte], color = colorstd[encontrado], linestyle='dotted') grafica.scatter(g_este,g_norte, color = 'blue', label = cadapunto) grafica.annotate(cadapunto,(g_este,g_norte), color='blue') grafica.set_xlabel('UTM_este') grafica.set_ylabel('UTM_norte') grafica.grid() grafica.set_title('Puntos Ubicados') plt.show()
-
5.3 Localiza por Trilateración – Sectores e Intervalos
Al incorporar en los algoritmos el uso de sectores, se obtienen resultados con menores errores de localización, estimados por trilateración y respecto a las coordenadas del gps.
Se realiza un enfoque para el área de vegetación, 'FIEC', para comparar resultados
Los resultados numéricos resumen obtenidos son:
puntos con medidas: 21 localizados con algoritmo: 20 , 95.2% Errores Promedio: +---------------+------+--------+--------+------+--------+--------+ | error+i*std | cant | % | trilat | gps | gps_dx | gps_dy | +---------------+------+--------+--------+------+--------+--------+ | interv_dentro | 19 | | | | | | | error_0 | 8 | 42.1% | 53.5 | 34.7 | 27.9 | 17.5 | | error_1 | 8 | 42.1% | 55.8 | 46.4 | 37.7 | 22.9 | | error_2 | 3 | 15.8% | 72.2 | 51.6 | 30.4 | 34.0 | | interv_fuera | 1 | | | | | | | error_0 | 0 | 0.0% | 0.0 | 0.0 | 0.0 | 0.0 | | error_1 | 1 | 100.0% | 24.1 | 27.0 | 27.0 | 1.6 | | error_2 | 0 | 0.0% | 0.0 | 0.0 | 0.0 | 0.0 | +---------------+------+--------+--------+------+--------+--------+
Se observa que se tienen casi todos los puntos del área de vegetación ubicados. 'FIEC130' no se incorpora a la lista, pues no disponde de valores de Rssi para 'gtwFCNM'.
El detalle de los errores de localización de cada punto son:
Errores localizacion Errores estimado: Cota-Trilatera-polígono y Trilatera_vs_GPS +---------+-------+-------+--------+-------+--------+--------+---------+---------+---------+ | punto | i*std | fuera | trilat | gps | gps_dx | gps_dy | sect_d1 | sect_d2 | sect_d3 | +---------+-------+-------+--------+-------+--------+--------+---------+---------+---------+ | FIEC101 | 1 | d2 | 24.13 | 27.01 | -26.96 | -1.61 | 0 | 0 | 0 | | FIEC102 | 1 | | 58.51 | 39.76 | 19.4 | -34.71 | 0 | 0 | 1 | | FIEC103 | 1 | | 30.16 | 59.15 | 50.75 | -30.38 | 0 | 0 | 1 | | FIEC104 | 2 | | 85.75 | 71.3 | -14.25 | -69.86 | 0 | 0 | 0 | | FIEC105 | 0 | | 2.92 | 22.87 | -22.11 | 5.83 | 0 | 0 | 0 | | FIEC106 | 0 | | 59.12 | 23.54 | -14.24 | -18.74 | 0 | 0 | 0 | | FIEC107 | 1 | | 26.94 | 15.29 | -5.65 | 14.21 | 0 | 0 | 0 | | FIEC108 | 1 | | 70.86 | 30.08 | -22.52 | 19.94 | 0 | 0 | 1 | | FIEC109 | 1 | | 61.94 | 21.15 | 16.68 | 13.0 | 0 | 0 | 1 | | FIEC110 | 0 | | 40.85 | 26.49 | -20.07 | 17.29 | 0 | 0 | 0 | | FIEC111 | 0 | | 72.84 | 17.7 | 1.5 | -17.63 | 0 | 0 | 1 | | FIEC112 | 0 | | 51.21 | 36.32 | 34.6 | 11.06 | 0 | 0 | 1 | | FIEC115 | 2 | | 78.12 | 10.24 | -9.45 | 3.96 | 1 | 0 | 0 | | FIEC116 | 1 | | 57.3 | 91.81 | 90.07 | 17.76 | 1 | 0 | 1 | | FIEC117 | 2 | | 52.61 | 73.24 | 67.54 | 28.33 | 1 | 0 | 1 | | FIEC120 | 0 | | 120.27 | 40.98 | -40.4 | -6.86 | 1 | 0 | 0 | | FIEC121 | 0 | | 49.52 | 50.94 | 40.54 | 30.84 | 1 | 0 | 0 | | FIEC122 | 1 | | 66.97 | 56.82 | 41.11 | 39.23 | 1 | 0 | 1 | | FIEC123 | 0 | | 30.95 | 58.98 | 49.64 | 31.85 | 1 | 0 | 0 | | FIEC124 | 1 | | 73.54 | 57.4 | 55.58 | -14.33 | 1 | 0 | 0 | +---------+-------+-------+--------+-------+--------+--------+---------+---------+---------+
En el detalle se observa que solo para 'FIEC101' se usó un intervalo extendido de la ecuación hacia 'd2' que es 'gtwFIEC' estimando que es por encontrarse en los límites de la sección de vegetación.
-
4.6 Rssi vs Distancia. Linealiza Sectores - Algoritmo Python
El efecto "sombra" añade el parámetro sector a la ecuación que se incorpora al resultado de la linealización en cada baliza.
La ecuación de una baliza se compone entonces de dos parámetros de selección:
sectoreintervalo.
Como ilustración se muestra la figura que tiene tres partes o ecuaciones:
ecuacion['s0']['r1']
ecuacion['s0']['r2']
ecuacion['s1']['r0']El primer parámetro para seleccionar la ecuación es el sector, 's0' y 's1', que en caso que sea un solo círculo se identifica como 's0'.
Dentro de cada sector, se mantiene el concepto de intervalos de distancia o radio. Se mantiene el concepto de la sección anterior, donde 'r0' corresponde a la linealización de todos los puntos en el sector. Cuando exiten sub-intervalos se usa 'r1', 'r2', etc para cada intervalo.
El número de ecuaciones corresponde al número de balizas y sectores establecidos para el análisis.
Se realizan cambios menores a la función pares_usar() de la librería girni para incorporar el parámetro sector, que al ser vacío
''funciona como fué descrito en las secciones anteriores.pares_usar(tabla,baliza, analiza,unabaliza, unsector ='', medida = 'rssi', modo = 'rx')
También se actualizaron los nombres de los archivos de entrada y salida para diferenciar de los resultados anteriores y disponer de los archivos para comparar con los resultados del método que solo usa intervalos.
Algoritmo en Python
# LoRa-Multipunto, Rssi vs distancia con mínimos cuadrados # linealización Rssi vs log10(distancia) # por Sectores e intervalos , Graficas '2D' # Girni 2020-10-07 propuesta: edelros@espol.edu.ec import numpy as np import pandas as pd import json import matplotlib.pyplot as plt import girni_lora_libreria as girni # INGRESO # archivos de entrada modo = 'rx' medida = 'rssi' arch_medidaubica = 'rsmP06_'+medida+'Ubica01sector1.txt' # archivos de salida arch_ecuaciones = 'rsmP07_ecuacionSector01.json' arch_medUbAtrib = 'rsmP07_'+medida+'UbicaUsarSector01.txt' # Analizar por segmentos analiza = {'gtwRECT':{'analizar' : 1, 'sector_ref' : ['FIEC112','FCNM110'], 's0':{'atipico_std' : 1, 'frontera' : [300], 'atipInterv_std': [2,2], 'p_amplia': 2, 'grp' : ['RECT','FIEC'], 'tip' : ['punto'], 'LOS' : [1] }, 's1':{'atipico_std' : 1, 'frontera' : [], 'atipInterv_std': [2], 'p_amplia': 2, 'grp' : ['FIEC','FCNM'], 'tip' : ['punto'], 'LOS' : [1,0] } }, 'gtwFIEC':{'analizar' : 1, 'sector_ref' : [], 's0':{'atipico_std' : 1, 'frontera' : [190], 'atipInterv_std': [1,1], 'p_amplia': 4, 'grp' : ['FIEC','FCNM'], 'tip' : ['punto'], 'LOS' : [0,1] } }, 'gtwFCNM':{'analizar' : 1, 'sector_ref' : [], 's0':{'atipico_std' : 1, 'frontera' : [235.0], 'atipInterv_std': [2,2], 'p_amplia': 4, 'grp' : ['FIEC','FCNM'], 'tip' : ['punto'], 'LOS' : [1,0] } } } baliza = {'d1':'gtwRECT', 'd2':'gtwFIEC', 'd3':'gtwFCNM'} # Parámetros de grafica tipograf = '2D' # '2D','' sin grafica escala = 'log' # 'normal','log' escalabase = 10 # 10 casicero = 1e-4 precision = 2 intersectar = 1 # 0:Falso, 1: Verdadero # PROCEDIMIENTO # leer datos tabla = pd.read_csv(arch_medidaubica, index_col='etiqueta') tabla = pd.DataFrame(tabla) baliza_key = list(baliza.keys()) baliza_val = list(baliza.values()) # Resultados de análisis ecuacion = {} eq_graf = {} # analiza datos hacia una baliza for unabaliza in analiza: donde = baliza_val.index(unabaliza) cualbaliza = baliza_key[donde] # Parámetros analizar = analiza[unabaliza]['analizar'] if analizar: # Crea ecuacion por baliza ecuacion[unabaliza] = {'sector_rad':[]} eq_graf[unabaliza] = {} # sectores sectores = [] sector_ref = analiza[unabaliza]['sector_ref'] tabla['sector_'+cualbaliza] = 0 # 'todos' predeterminado # coordenadas baliza para identificar ángulo b_este = tabla['c_este'][unabaliza] b_norte = tabla['c_norte'][unabaliza] # sectores por puntos de referencia if len(sector_ref)>0: for cadapunto in sector_ref: p_este = tabla['c_este'][cadapunto] p_norte = tabla['c_norte'][cadapunto] dx = p_este - b_este dy = p_norte - b_norte theta = np.arctan2(dy,dx) if theta<0 and dx<0: theta = theta + 2*np.pi sectores.append(theta) sectores = np.array(sectores) ordenar = np.argsort(sectores) sectores = list(sectores[ordenar]) nsectores = len(sectores) # clasifica puntos por sector particular for cadapunto in tabla.index: dentrosector = 0 p_este = tabla['c_este'][cadapunto] p_norte = tabla['c_norte'][cadapunto] dx = p_este-b_este dy = p_norte-b_norte theta = np.arctan2(dy,dx) if theta<0 and dx<0: theta = theta + 2*np.pi for j in range(0,nsectores-1,1): if theta>sectores[j] and theta<sectores[j+1]: dentrosector = j+1 tabla.loc[cadapunto,'sector_'+cualbaliza] = dentrosector # ecuacion por baliza y sector ecuacion[unabaliza]['sector_rad'] = sectores nsectores = len(sectores) if nsectores == 0: nsectores = 1 for cadasector in range(0,nsectores,1): unsector = 's'+str(cadasector) ecuacion[unabaliza][unsector] = {} eq_graf[unabaliza][unsector] = {} # ecuación con todos los puntos como referencia [pares,par_etiqueta] = girni.pares_usar(tabla,baliza,analiza, unabaliza,unsector, medida,modo) xi = pares[:,0] yi = pares[:,1] n_xi = len(xi) # coeficiente de correlación correlacion = np.corrcoef(xi,yi)[0,1] # minimos cuadrados ecuacion0 = girni.linealiza_lstsq(xi,yi) # selecciona atipicos de todos los puntos atipico_std = analiza[unabaliza][unsector]['atipico_std'] alpha = ecuacion0['alpha'] beta = ecuacion0['beta'] fdist0 = lambda d: -10*alpha*(np.log10(d))+beta yi0 = fdist0(xi) dyi0std = ecuacion0['error_std'] dyi0 = yi - yi0 atipicos = np.abs(dyi0) >= dyi0std*atipico_std xi0_e = xi[atipicos] yi0_e = yi[atipicos] etiq0_e = par_etiqueta[atipicos] unintervalo = 'r0' # todos # para exportar hacia archivo o gráfica ecuacion[unabaliza][unsector] = {unintervalo: ecuacion0 } ecuacion[unabaliza][unsector][unintervalo]['correlacion'] = correlacion eq_graf[unabaliza][unsector] = {unintervalo: {'xi_graf' : xi, 'yi_graf' : yi, 'etiqueta' : par_etiqueta, 'linea' : yi0, 'atipicos' : [xi0_e,yi0_e], 'atip_etiq': etiq0_e} } # Intervalos radiales en sector intervalo = [np.min(xi),np.max(xi)] frontera = analiza[unabaliza][unsector]['frontera'] if len(frontera)>0: # revisar si frontera esta dentro intervalo frontera = np.array(frontera, dtype=float) revisar = (frontera>=np.min(xi)) & (frontera<=np.max(xi)) enintervalo = list(frontera[revisar]) intervalo.extend(enintervalo) intervalo = np.array(intervalo) ordenar = np.argsort(intervalo) intervalo = intervalo[ordenar] n_intervalo = len(intervalo) # analizar cada subintervalo p_inicio = 0 p_desde = 0 p_amplia = analiza[unabaliza][unsector]['p_amplia'] atipIntv_std = analiza[unabaliza][unsector]['atipInterv_std'] for i in range(0,n_intervalo-1,1): i_eq = 'r' + str(i+1) # puntos en subintervalo [a,b] a = intervalo[i] b = intervalo[i+1] subintervalo = (xi >= a) & (xi <= b) xi_sub = xi[subintervalo] yi_sub = yi[subintervalo] n_sub = len(xi_sub) etiq_sub = par_etiqueta[p_inicio:p_inicio + n_sub] # amplia sub-intervalo, mejora intersecta rectas detras = p_inicio retrocede = detras if detras > p_amplia: retrocede = p_amplia delante = n_xi - (p_inicio+n_sub)# -1) avanza = delante if delante >= p_amplia: avanza = p_amplia p_desde = p_inicio - retrocede p_hasta = p_inicio + (n_sub) + avanza p_inicio = p_inicio + (n_sub-1) # subintervalo, amplia puntos xi_a = xi[p_desde:p_hasta] yi_a = yi[p_desde:p_hasta] etiq_a = par_etiqueta[p_desde:p_hasta] # coeficiente de correlación correlacion1 = np.corrcoef(xi_a,yi_a)[0,1] # analiza subintervalo ecuacion1 = girni.linealiza_lstsq(xi_a,yi_a) ecuacion[unabaliza][unsector][i_eq] = ecuacion1 ecuacion[unabaliza][unsector][i_eq]['correlacion'] = correlacion1 # atipicos del subintervalo extendido alpha = ecuacion1['alpha'] beta = ecuacion1['beta'] fdist1 = lambda d: -10*alpha*(np.log10(d))+beta yi1 = fdist1(xi_a) dyi1std = ecuacion1['error_std'] atipicos = np.zeros(len(xi_a),dtype=bool) atipico_std = analiza[unabaliza][unsector]['atipInterv_std'][i] dyi1 = yi_a - yi1 if np.abs(dyi1std) > casicero: atipicos = np.abs(dyi1) >= dyi1std*atipico_std xi1_e = xi_a[atipicos] yi1_e = yi_a[atipicos] etiq1_e = etiq_a[atipicos] # para gráfica, atipicos sin extender puntos atipicos_sub = (xi1_e >= a) & (xi1_e<=b) xi_sub1_e = xi1_e[atipicos_sub] yi_sub1_e = yi1_e[atipicos_sub] etiq_sub1e = etiq1_e[atipicos_sub] eq_graf[unabaliza][unsector][i_eq] = {'atipicos': [xi_sub1_e,yi_sub1_e], 'atip_etiq':etiq_sub1e} # subintervalo sin atipicos if len(xi1_e)>0: atipicoNo = np.abs(dyi1) <= dyi1std*atipico_std xi2 = xi_a[atipicoNo] yi2 = yi_a[atipicoNo] etiq2 = etiq_a[atipicoNo] # coeficiente de correlación correlacion2 = np.corrcoef(xi2,yi2)[0,1] ecuacion2 = girni.linealiza_lstsq(xi2,yi2) # actualiza ecuación sin atipicos intervaloy intervalox = ecuacion1['intervalox'] ecuacion2['intervalox'] = intervalox.copy() alpha = ecuacion2['alpha'] beta = ecuacion2['beta'] fdist = lambda d: -10*alpha*(np.log10(d))+beta intervaloy = fdist(intervalox) ordenar = np.argsort(intervaloy) intervaloy = list(intervaloy[ordenar]) ecuacion2['intervaloy'] = intervaloy ecuacion[unabaliza][unsector][i_eq] = ecuacion2 ecuacion[unabaliza][unsector][i_eq]['correlacion'] = correlacion2 # Revisar frontera entre subintervalos, # calcula intersección de rectas como nueva frontera interv_calc = np.copy(intervalo) if len(intervalo) >2 and intersectar==1 : for i in range(0,n_intervalo-2,1): ai = 'r' + str(i+1) bi = 'r' + str(i+2) ma = ecuacion[unabaliza][unsector][ai]['alpha'] ba = ecuacion[unabaliza][unsector][ai]['beta'] mb = ecuacion[unabaliza][unsector][bi]['alpha'] bb = ecuacion[unabaliza][unsector][bi]['beta'] # punto de intersección o cruce cruzanx = 10**((bb-ba)/(10*(mb-ma))) dfrontera = frontera - cruzanx # cruce dentro de intervalo de ecuacion if cruzanx > intervalo[-1]: cruzanx = intervalo[-1] if cruzanx < intervalo[0]: cruzanx = intervalo[0] interv_calc[i+1] = cruzanx # para grafica evalua cada subintervalo sin atipicos n_interv_calc = len(interv_calc) for i in range(0,n_interv_calc-1,1): i_eq = 'r'+str(i+1) a = interv_calc[i] b = interv_calc[i+1] subintervalo = (xi >= a) & (xi <= b) xi_sub = xi[subintervalo] yi_sub = yi[subintervalo] xi_graf = np.copy(xi[subintervalo]) if not(a in xi_graf): xi_graf = np.concatenate(([a],xi_graf),axis=0) if not(b in xi_graf): xi_graf = np.concatenate((xi_graf,[b]),axis=0) # Evalua subintervalo con la ecuacion sin atipicos alpha = ecuacion[unabaliza][unsector][i_eq]['alpha'] beta = ecuacion[unabaliza][unsector][i_eq]['beta'] fdist = lambda d: -10*alpha*(np.log10(d))+beta yi1_sub = fdist(xi_sub) yi_graf = fdist(xi_graf) eq_graf[unabaliza][unsector][i_eq]['xi_graf'] = xi_graf eq_graf[unabaliza][unsector][i_eq]['yi_graf'] = yi_graf a = np.round(np.min([xi_graf]),precision) b = np.round(np.max([xi_graf]),precision) ecuacion[unabaliza][unsector][i_eq]['intervalox'] = [a,b] ay = np.round(np.min([yi_graf]),precision) by = np.round(np.max([yi_graf]),precision) ecuacion[unabaliza][unsector][i_eq]['intervaloy'] = [ay,by] # SALIDA for unabaliza in ecuacion: for unsector in ecuacion[unabaliza]: if unsector == 'sector_rad': print('baliza: ',unabaliza) print(' sectores radianes: ',ecuacion[unabaliza]['sector_rad']) if unsector != 'sector_rad': for i_eq in ecuacion[unabaliza][unsector]: unintervalo = ecuacion[unabaliza][unsector][i_eq]['intervalox'] unintervaloy = ecuacion[unabaliza][unsector][i_eq]['intervaloy'] error_medio = ecuacion[unabaliza][unsector][i_eq]['error_medio'] error_std = ecuacion[unabaliza][unsector][i_eq]['error_std'] eq_latex = ecuacion[unabaliza][unsector][i_eq]['eq_latex'] errorx_medio = ecuacion[unabaliza][unsector][i_eq]['errorx_medio'] errorx_std = ecuacion[unabaliza][unsector][i_eq]['errorx_std'] correlacion = ecuacion[unabaliza][unsector][i_eq]['correlacion'] print(' [sector][intervalo]: ',unsector,',',i_eq) print(' ' + eq_latex) print(' ','intervalox: ',np.round(unintervalo,precision)) print(' ','intervaloy: ',np.round(unintervaloy,precision)) print(' correlación: ',np.round(correlacion,precision)) print(' |error_rssi| promedio: ',np.round(error_medio,precision), ' , std:',np.round(error_std,precision)) print(' |error_dist| promedio: ',np.round(errorx_medio,precision), ' , std:',np.round(errorx_std,precision)) print() # salida hacia archivo with open(arch_ecuaciones, 'w') as outfile: json.dump(ecuacion, outfile) tabla.to_csv(arch_medUbAtrib) # GRAFICA # Referencias para gráfica grupo = ['FIEC' ,'FCNM' ,'RECT','CIRC'] colores = ['green','orange','grey','magenta'] tipo = ['punto','1m' ,'gtw','dispositivo'] marcas = [ 'o','D' ,'D' ,'*' ] if tipograf=='2D': for unabaliza in ecuacion: for unsector in ecuacion[unabaliza]: if not(unsector=='sector_rad'): figura,grafica = plt.subplots() if escala == 'log': grafica.set_xscale(escala,base=escalabase) # todos los puntos unintervalo = 'r0' xi = eq_graf[unabaliza][unsector][unintervalo]['xi_graf'] yi = eq_graf[unabaliza][unsector][unintervalo]['yi_graf'] etiqueta = eq_graf[unabaliza][unsector][unintervalo]['etiqueta'] grafica.scatter(xi,yi,marker='.') m = len(xi) for i in range(0,m,1): grafica.annotate(etiqueta[i],(xi[i],yi[i])) # linealizado con todos los puntos fdtxt = ecuacion[unabaliza][unsector][unintervalo]['eq_latex'] yi0 = eq_graf[unabaliza][unsector][unintervalo]['linea'] a = np.round(np.min([xi]),precision) b = np.round(np.max([xi]),precision) eq_texto = fdtxt +' ; ['+ str(a) +','+ str(b)+']' grafica.plot(xi,yi0,label=eq_texto,linestyle='dotted') # linealizado por subintervalo eq_interv = list(ecuacion[unabaliza][unsector].keys()) eq_interv.pop(0) n_intervalo = len(eq_interv) for i_eq in eq_interv: fdtxt = ecuacion[unabaliza][unsector][i_eq]['eq_latex'] grtxt = ecuacion[unabaliza][unsector][i_eq]['eqg_latex'] xi_graf = eq_graf[unabaliza][unsector][i_eq]['xi_graf'] yi_graf = eq_graf[unabaliza][unsector][i_eq]['yi_graf'] a = np.round(np.min([xi_graf]),precision) b = np.round(np.max([xi_graf]),precision) eq_texto = fdtxt+' ; ['+str(a)+','+str(b)+']' grafica.plot(xi_graf,yi_graf, label=eq_texto) # atipicos marcados en subintervalo [xi1_e,yi1_e] = eq_graf[unabaliza][unsector][i_eq]['atipicos'] etiq1_e = eq_graf[unabaliza][unsector][i_eq]['atip_etiq'] grafica.scatter(xi1_e,yi1_e, color='red') m = len(etiq1_e) for i in range(0,m,1): grafica.annotate(etiq1_e[i], (xi1_e[i],yi1_e[i]), color='red') # lineas de frontera grafica.axvline(a, color='lightblue') valor_frontera = str(np.round(a,precision)) grafica.annotate(valor_frontera, (a,np.max([yi,yi0])), color='lightblue') grafica.axvline(b, color='lightblue') valor_frontera = str(np.round(b,precision)) grafica.annotate(valor_frontera, (b,np.max([yi,yi0])), color='lightblue') # etiquetas y títulos grafica.legend() grafica.set_ylabel(medida+'_'+modo) grafica.set_xlabel('distancia') grafica.grid(True,linestyle='dotted', axis='x', which='both') untitulo = unabaliza+'_'+unsector+': ' untitulo = untitulo+medida+'_'+modo+' vs distancia' grafica.set_title(untitulo) plt.show()
-
4.5 Rssi vs Distancia. Linealiza Sectores
Para la baliza de Rectorado, hacia el área de estudio se encuentra el edificio Biblioteca que es suficientemente grande para bloquear parcialmente la señal y genera un efecto "sombra".
Biblioteca ESPOL En la vista superior del área de mediciones, a la derecha se ubica el rectorado en la cima del cerro, la Biblioteca en una parte intermedia, y el área de vegetacion FIEC se ubica en la parte baja del cerro hacia la izquierda de la imagen.
La "sombra" en la parte de vegetación se delimita en el sector formado por el ángulo formado entre las líneas de baliza y punto de referencia.
En la imagen se identifica como sector 's0' al círculo completo, como si no hubiese sectorización, para luego realizar cortes para el sector 's1' de la "sombra".
El sector 's1', del ejemplo, usa los puntos de referencia de inicio y fin en el sentido de las manecillas del reloj:
-
'FIEC112'usado como inicio, pues la línea entre Rectorado y el punto FIEC112 pasa por el borde del edificio de Biblioteca.
-'FCNM110'es usado punto final del sector al delimitar el área de vegetación bajo estudio entre FIEC y FCNM desde RECTorado.En cada sector, se usan también intervalos por distancia desde la baliza descrita en la sección anterior.
El asunto con sectores se presenta principalmente con la baliza 'gtwRECT', siendo el caso referencia para el desarrollo del algoritmo. El sector añade una variable y nivel al diccionario de la ecuación.
baliza: ‘gtwRECT’ , Sector: 's0'
Para el primer sector 's0' sin "somlbra" del edificion biblioteca, también se aplica un subintervalo con frontera buscada desde los 300 m y determinada con las fórmulas en un radio de 359 m.
baliza: ‘gtwFIEC’ , Sector: 's1'
Para el siguiente sector 's1' se observa el efecto de "sombra" del edificio Biblioteca, se presenta una recuperación de señal al alejar de la baliza de RECTorado y el edificio Biblioteca.
En la gráfica, la baliza de RECTorado se encuentra a la izquierda, pues el eje distancia aumenta de valor hacia la derecha.
Parámetros de sector
La variable
'sector_ref'contiene los puntos de referencia para un nuevo sector 's1'. Si 'sector_ref' no contiene valores[]se asume que una formula aplica todo el círculo predeterminada como sector 's0'.Los puntos de referencia se convierten a su equivalente en radianes usando las coordenadas utm este y norte, valores en radianes que pasan a formar parte de diccionario de ecuaciones como
'sector_rad'usados luego para evaluar la ecuación.# Analizar por segmentos analiza = {'gtwRECT':{'analizar' : 1, 'sector_ref' : ['FIEC112','FCNM110'], 's0':{'atipico_std' : 1, 'frontera' : [300], 'atipInterv_std': [2,2], 'p_amplia': 2, 'grp' : ['RECT','FIEC'], 'tip' : ['punto'], 'LOS' : [1] }, 's1':{'atipico_std' : 1, 'frontera' : [], 'atipInterv_std': [2], 'p_amplia': 2, 'grp' : ['FIEC','FCNM'], 'tip' : ['punto'], 'LOS' : [1,0] } }, 'gtwFIEC':{'analizar' : 1, 'sector_ref' : [], 's0':{'atipico_std' : 1, 'frontera' : [190], 'atipInterv_std': [1,1], 'p_amplia': 4, 'grp' : ['FIEC','FCNM'], 'tip' : ['punto'], 'LOS' : [0,1] } }, 'gtwFCNM':{'analizar' : 1, 'sector_ref' : [], 's0':{'atipico_std' : 1, 'frontera' : [235.0], 'atipInterv_std': [2,2], 'p_amplia': 4, 'grp' : ['FIEC','FCNM'], 'tip' : ['punto'], 'LOS' : [1,0] } } }
Resultados en otras balizas
El caso de la baliza 'gtwFIEC' no tiene observaciones de "sombra" para el área de vegetación que sean de tamaño considerable, por lo que se mantiene solo la división por segmentos.
Los resultados no varían gtwFCNM, pues no se ha aplicado este concepto de "sombra" para estas balizas, tan solo la división por subintevalos causadas por cambios de entorno.
Los resultados del algoritmo para usar en localización :
baliza: gtwRECT sectores radianes: [2.6552133371486635, 2.7999991438538827] [sector][intervalo]: s0 , r0 $ rssi = -10(5.148)log_{10}(d)+(15.043) $ intervalox: [138.15 482.74] intervaloy: [-123.1 -93.18] correlación: -0.93 |error_rssi| promedio: 2.88 , std: 3.55 |error_dist| promedio: 39.64 , std: 44.91 [sector][intervalo]: s0 , r1 $ rssi = -10(5.001)log_{10}(d)+(11.842) $ intervalox: [138.15 359.35] intervaloy: [-115.97 -95.2 ] correlación: -0.9 |error_rssi| promedio: 4.18 , std: 4.75 |error_dist| promedio: 42.76 , std: 44.83 [sector][intervalo]: s0 , r2 $ rssi = -10(7.191)log_{10}(d)+(67.794) $ intervalox: [359.35 482.74] intervaloy: [-125.18 -115.97] correlación: -0.94 |error_rssi| promedio: 2.11 , std: 2.3 |error_dist| promedio: 25.23 , std: 28.25 [sector][intervalo]: s1 , r0 $ rssi = -10(-8.383)log_{10}(d)+(-337.543) $ intervalox: [374.66 443.48] intervaloy: [-124.2 -113.71] correlación: 0.66 |error_rssi| promedio: 1.79 , std: 2.2 |error_dist| promedio: 20.48 , std: 25.57 [sector][intervalo]: s1 , r1 $ rssi = -10(-8.383)log_{10}(d)+(-337.543) $ intervalox: [374.66 443.48] intervaloy: [-121.79 -115.65] correlación: 0.66 |error_rssi| promedio: 1.79 , std: 2.2 |error_dist| promedio: 20.48 , std: 25.57 baliza: gtwFIEC sectores radianes: [] [sector][intervalo]: s0 , r0 $ rssi = -10(4.908)log_{10}(d)+(1.406) $ intervalox: [ 52.54 397.15] intervaloy: [-129.13 -86.99] correlación: -0.9 |error_rssi| promedio: 4.84 , std: 5.56 |error_dist| promedio: 41.43 , std: 54.09 [sector][intervalo]: s0 , r1 $ rssi = -10(4.263)log_{10}(d)+(-11.269) $ intervalox: [ 52.54 166.14] intervaloy: [-105.93 -84.62] correlación: -0.9 |error_rssi| promedio: 2.92 , std: 3.33 |error_dist| promedio: 20.19 , std: 24.67 [sector][intervalo]: s0 , r2 $ rssi = -10(6.09)log_{10}(d)+(29.295) $ intervalox: [166.14 397.15] intervaloy: [-128.98 -105.93] correlación: -0.92 |error_rssi| promedio: 2.71 , std: 3.17 |error_dist| promedio: 25.03 , std: 29.93 baliza: gtwFCNM sectores radianes: [] [sector][intervalo]: s0 , r0 $ rssi = -10(5.403)log_{10}(d)+(8.423) $ intervalox: [ 27.74 364.71] intervaloy: [-132.09 -81.62] correlación: -0.93 |error_rssi| promedio: 4.59 , std: 5.48 |error_dist| promedio: 41.33 , std: 48.24 [sector][intervalo]: s0 , r1 $ rssi = -10(4.795)log_{10}(d)+(-3.65) $ intervalox: [ 27.74 238.81] intervaloy: [-117.68 -72.85] correlación: -0.91 |error_rssi| promedio: 5.96 , std: 6.51 |error_dist| promedio: 52.15 , std: 64.38 [sector][intervalo]: s0 , r2 $ rssi = -10(9.027)log_{10}(d)+(96.989) $ intervalox: [238.81 364.71] intervaloy: [-134.28 -117.68] correlación: -0.87 |error_rssi| promedio: 2.81 , std: 3.43 |error_dist| promedio: 20.26 , std: 25.06 >>> -
5.2 Localiza por Trilateración – Intervalos. Algoritmo Python
Algoritmo en Python
# Localización por trilateración # Ubicación de puntos en el mapa # Girni 2020-10-07 propuesta: edelros@espol.edu.ec import numpy as np import sympy as sym import matplotlib.pyplot as plt import pandas as pd import json import girni_lora_libreria as girni # INGRESO # revisar parametros al inicio modo = 'rx' medida = 'rssi' # archivos de entrada arch_medUbAtrib = 'rsmP06_'+medida+'UbicaUsar01.txt' arch_ecuaciones = 'rsmP06_ecuaciones01.json' # archivos de salida arch_trilatera = 'rsmP07_'+medida+'trilatera01.txt' # grupo a procesar mostrargrp = ['FIEC'] #['FIEC','FCNM','RECT','CIRC'] mostrartip = ['punto'] mostrarfuera = [0,1] #[0,1] tolera_error = 2 # decimales en error # Referencias baliza = {'d1':'gtwRECT', 'd2':'gtwFIEC', 'd3':'gtwFCNM'} # PROCEDIMIENTO baliza_key = list(baliza.keys()) baliza_val = list(baliza.values()) # leer datos with open(arch_ecuaciones) as json_file: ecuacion = json.load(json_file) tabla = pd.read_csv(arch_medUbAtrib, index_col='etiqueta') tabla = pd.DataFrame(tabla) # coordenadas de baliza coordbaliza = {} for cadabaliza in baliza: cualbaliza = baliza[cadabaliza] coordbaliza[cadabaliza] = [tabla['c_este'][cualbaliza]] coordbaliza[cadabaliza].append(tabla['c_norte'][cualbaliza]) # inicializa ubicados, encontrados, intervalos fuera, errores # ubicados tiene [distancia, distancia+1s, distancia +2s,total] ubicados = np.zeros(shape=(3,4),dtype=int) tabla['interv_fuera'] = '' tabla['encontrado'] = np.nan # suma de errores sumatrilat = np.zeros(shape=(2,4),dtype=float) sumagps = np.zeros(shape=(2,4),dtype=float) sumagpsdx = np.zeros(shape=(2,4),dtype=float) sumagpsdy = np.zeros(shape=(2,4),dtype=float) # distancias y errores a cada baliza for unabaliza in baliza_val: donde = baliza_val.index(unabaliza) cualbaliza = baliza_key[donde] columna = medida+'_'+modo+'_'+cualbaliza tabla['dist0_'+cualbaliza] = np.nan tabla['dist1_'+cualbaliza] = np.nan tabla['dist2_'+cualbaliza] = np.nan # distancia estimada por punto for cadapunto in tabla.index: p_rssi = tabla[columna][cadapunto] esgrupo = tabla['grupo'][cadapunto] estipo = tabla['tipo'][cadapunto] # ecuacion a usar ecuacion_rssi = ecuacion[unabaliza] if not(np.isnan(p_rssi)): if (esgrupo in mostrargrp) and (estipo in mostrartip): ubicados[2,donde] = ubicados[2,donde] + 1 dist = girni.dist_rssi(p_rssi,ecuacion_rssi) [distancia,e_mean,e_1std,e_2std,interv_fuera] = dist dist1 = distancia + e_1std dist2 = distancia + e_2std tabla.loc[cadapunto,'dist0_'+cualbaliza] = distancia tabla.loc[cadapunto,'dist1_'+cualbaliza] = dist1 tabla.loc[cadapunto,'dist2_'+cualbaliza] = dist2 if interv_fuera>0: if len(tabla['interv_fuera'][cadapunto])==0: separador='' else: separador=',' tabla.loc[cadapunto,'interv_fuera'] = tabla['interv_fuera'][cadapunto]+separador+cualbaliza procesados = np.max(ubicados[2,:]) # trilateración. estima localización cada punto # inicializa cada punto a i*std for i in range(0,3,1): tabla['trilat_este_'+str(i)] = np.nan tabla['trilat_norte_'+str(i)] = np.nan tabla['trilat_error_'+str(i)] = np.nan tabla['ubicado_'+str(i)] = np.nan tabla['trilat_gps_error_'+str(i)] = np.nan tabla['trilat_gps_dx_'+str(i)] = np.nan tabla['trilat_gps_dy_'+str(i)] = np.nan # procesa cada punto for cadapunto in tabla.index: esgrupo = tabla['grupo'][cadapunto] if esgrupo in mostrargrp: # coordenadas GPS del punto c_este = tabla['c_este'][cadapunto] c_norte = tabla['c_norte'][cadapunto] # cada punto a i*std punto = {} for i in range(0,3,1): punto[i] = {} cuenta = 0 for unabaliza in baliza: punto[i][unabaliza] = tabla['dist'+str(i)+'_'+unabaliza][cadapunto] # si existe distancia if not(np.isnan(punto[i][unabaliza])): cuenta = cuenta + 1 # Hay 3 mediciones a baliza if cuenta==3: localiza = girni.trilatera(punto[i], coordbaliza, tolera = 10e-4) baricentro = localiza['baricentro'] barerror = localiza['barerror'] poligono = localiza['poligono'] sumabaricentro = np.sum(baricentro) # Hay coordenadas de baricentro if not(np.isnan(sumabaricentro)): tabla.loc[cadapunto,'trilat_este_'+str(i)] = baricentro[0] tabla.loc[cadapunto,'trilat_norte_'+str(i)] = baricentro[1] tabla.loc[cadapunto,'trilat_error_'+str(i)] = np.round(barerror,tolera_error) # error trilatera hacia gps dx = baricentro[0] - c_este dy = baricentro[1] - c_norte error_gps = np.sqrt(dx**2+dy**2) tabla.loc[cadapunto,'trilat_gps_error_'+str(i)] = np.round(error_gps,tolera_error) tabla.loc[cadapunto,'trilat_gps_dx_'+str(i)] = np.round(dx,tolera_error) tabla.loc[cadapunto,'trilat_gps_dy_'+str(i)] = np.round(dy,tolera_error) # cuenta error si encuentra primera vez cond1 = np.isnan(tabla['encontrado'][cadapunto]) cond2 = len(tabla['interv_fuera'][cadapunto])>0 k = int(cond2) # intervalo extendido if cond1: tabla.loc[cadapunto,'encontrado'] = i ubicados[k,i] = ubicados[k,i]+1 sumatrilat[k,i] = sumatrilat[k,i] + np.round(barerror,tolera_error) sumagps[k,i] = sumagps[k,i] + np.round(error_gps,tolera_error) sumagpsdx[k,i] = sumagpsdx[k,i] + np.round(np.abs(dx),tolera_error) sumagpsdy[k,i] = sumagpsdy[k,i] + np.round(np.abs(dy),tolera_error) # contabiliza errores de localizados ubicasuma = np.zeros(2) ubicasumporc = np.zeros(2) ubicaporc = np.zeros(shape=(2,3)) for k in range(0,2,1): ubicados[k,3] = int(np.sum(ubicados[k,0:3])) ubicasuma[k] = np.sum(ubicados[k,0:3]) ubicasumporc[k] = np.round(100*ubicasuma[k]/ubicados[k,3],1) for i in range(0,3,1): ubicaporc[k,i] = np.round(100*ubicados[k,i]/ubicasuma[k],1) if ubicados[k,i]>0: sumatrilat[k,i] = np.round(sumatrilat[k,i]/ubicados[k,i],1) sumagps[k,i] = np.round(sumagps[k,i]/ubicados[k,i],1) sumagpsdx[k,i] = np.round(sumagpsdx[k,i]/ubicados[k,i],1) sumagpsdy[k,i] = np.round(sumagpsdy[k,i]/ubicados[k,i],1) # SALIDA print('Errores localizacion') import prettytable as ptt print('Errores estimado: Cota-Trilatera-polígono y Trilatera_vs_GPS') mostrar = ptt.PrettyTable(['punto','i*std','fuera','trilat', 'gps','gps_dx','gps_dy', 'u_d1','u_d2','u_d3']) for cadapunto in tabla.index: # selecciona puntos a mostrar encontrado = tabla['encontrado'][cadapunto] esgrupo = tabla['grupo'][cadapunto] interv_fuera = int(len(tabla['interv_fuera'][cadapunto])!=0) cond1 = not(np.isnan(encontrado)) cond2 = esgrupo in mostrargrp cond3 = interv_fuera in mostrarfuera if cond1 and cond2 and cond3: encontrado = int(encontrado) mostrar.add_row([cadapunto, str(int(tabla['encontrado'][cadapunto])), tabla['interv_fuera'][cadapunto], str(tabla['trilat_error_'+str(encontrado)][cadapunto]), str(tabla['trilat_gps_error_'+str(encontrado)][cadapunto]), str(tabla['trilat_gps_dx_'+str(encontrado)][cadapunto]), str(tabla['trilat_gps_dy_'+str(encontrado)][cadapunto]), tabla['usar_d1'][cadapunto], tabla['usar_d2'][cadapunto], tabla['usar_d3'][cadapunto]]) print(mostrar) print('puntos con medidas: ',procesados) loc_alg = np.round(100*np.sum(ubicados[:,3])/procesados,1) print('localizados con algoritmo: ', np.sum(ubicados[:,3]), ' , '+str(loc_alg)+'%' ) print('Errores Promedio:') resumen = ptt.PrettyTable(['error+i*std','cant','%', 'trilat','gps','gps_dx', 'gps_dy']) for k in range(0,2,1): if k == 0: texto = 'dentro' if k == 1: texto = 'fuera' resumen.add_row(['interv_'+texto,ubicados[k,3],'', '','','','']) for i in range(0,3,1): resumen.add_row(['error_'+str(i),ubicados[k,i], str(ubicaporc[k,i])+'%', sumatrilat[k,i],sumagps[k,i], sumagpsdx[k,i],sumagpsdy[k,i], ]) print(resumen) # salida hacia archivo tabla.to_csv(arch_trilatera) # Grafica ubicados ------------ # Referencias para gráfica grupo = ['FIEC' ,'FCNM' ,'RECT','CIRC'] colores = ['green','orange','grey','magenta'] tipo = ['punto','1m' ,'gtw','dispositivo'] marcas = [ 'o','D' ,'D' ,'*' ] colorstd = ['lightblue', 'lightgreen','orange'] colorlin = ['lightblue', 'lightgreen','orange'] figura,grafica = plt.subplots() # balizas for unabaliza in coordbaliza: g_este = coordbaliza[unabaliza][0] g_norte = coordbaliza[unabaliza][1] grafica.scatter(g_este,g_norte, color = 'red', marker = 'D', label = cadapunto) grafica.annotate(unabaliza, (g_este,g_norte)) # Puntos for cadapunto in tabla.index: g_este = tabla['c_este'][cadapunto] g_norte = tabla['c_norte'][cadapunto] # selecciona puntos a mostrar encontrado = tabla['encontrado'][cadapunto] esgrupo = tabla['grupo'][cadapunto] interv_fuera = int(len(tabla['interv_fuera'][cadapunto])!=0) cond1 = not(np.isnan(encontrado)) cond2 = esgrupo in mostrargrp cond3 = interv_fuera in mostrarfuera if cond1 and cond2 and cond3: encontrado = int(encontrado) p_este = tabla['trilat_este_'+str(encontrado)][cadapunto] p_norte = tabla['trilat_norte_'+str(encontrado)][cadapunto] grafica.scatter(p_este,p_norte, color = colorstd[encontrado], label = cadapunto) grafica.plot([p_este,g_este], [p_norte,g_norte], color = colorstd[encontrado], linestyle='dotted') grafica.scatter(g_este,g_norte, color = 'blue', label = cadapunto) grafica.annotate(cadapunto,(g_este,g_norte), color='blue') grafica.set_xlabel('UTM_este') grafica.set_ylabel('UTM_norte') grafica.grid() grafica.set_title('Puntos Ubicados') plt.show()
-
4.2 Rssi vs distancia. Linealiza - función Python
Para el procesamiento de los datos se incorpora una función a girni_lora_libreria para realizar la linealización por mínimos cuadrados.
linealiza_lstsq(xi,yi,digitos = 3)La función se encarga de convertir el eje x en log10(x), asi como construir las ecuaciones en la forma numérica lambda, latex y un diccionario con los parámetros de la ecuación.
Datos de ingreso
Los datos de ingreso son xi y yi para cada eje, la variable dígitos establece los decimales a usar en la expresión en formato latex.
Datos de Salida
El resultado es un diccionario con los intervalos de los valores obtenidos para cada eje, la pendiente de la recta, el |error| promedio, la desviación estándar error_std, la ecuación en formato latex.
unaecuacion = {'alpha' : alpha, 'beta' : beta, 'eq_latex': fdtxt0, 'intervalox' : [np.min(xi),np.max(xi)], 'error_medio': dyi0mean, 'error_std' : dyi0std, 'eqg_latex' : grtxt0, 'intervaloy' : [np.min(yi),np.max(yi)], 'errorx_medio': dxi0mean, 'errorx_std' : dxi0std, }El resultado se puede escribir en un archivo en formato json. La ecuación se recupera desde el archivo con lo que se puede volver a construir la función en la lambda para evaluación numérica.
Procedimiento
Se desarrolla principalmente usando la función numpy.linalg.lstsq.
La relación rssi vs distancia usa log10(xi), por lo que se incluye esta operación antes de aplicar mínimos cuadrados.
Obtenidos los parámetros, se da el formato de la expresión acorde al modelo básico de pérdidas en espacio libre en latex para mostrar como etiqueta en las gráficas.
RSSI(d) = -10 \alpha \log_{10} (d) + P_{0}La función se usa para cada baliza, y en varios segmentos, para observar los posibles resultados, también se incorporan los valores de errores,
Algoritmo en Python
def linealiza_lstsq(xi,yi,digitos = 3): ''' usa minimos cuadrados para entregar la ecuacion digitos: usados en expresion latex ''' unaecuacion = {} # Eje x en log10() xilog = np.log10(xi) n = len(xi) # mínimos cuadrados (least square), # distancia vs medida A = np.vstack([xilog, np.ones(n)]).T [m0, b0] = np.linalg.lstsq(A, yi, rcond=None)[0] alpha = -m0/10 beta = b0 # ecuaciones expresion rssi(d) fdist0 = lambda d: -10*alpha*(np.log10(d))+beta fdtxt0 = r'$ rssi = -10(' + str(np.round(alpha,digitos)) fdtxt0 = fdtxt0 + ')log_{10}(d)+(' fdtxt0 = fdtxt0 + str(np.round(beta,digitos))+') $' # Errores respecto a rssi(d) yi0 = fdist0(xi) dyi0 = yi - yi0 dyi0mean = np.mean(np.abs(dyi0)) dyi0std = np.std(dyi0, dtype=np.float64) # ecuaciones expresion d(rssi) grssi0 = lambda rssi: 10**((beta-rssi)/(10*alpha)) grtxt0 = r"$ d = 10^{(" + str(np.round(beta,digitos)) + ' - ' grtxt0 = grtxt0 + 'rssi)/' + '(10('+str(np.round(alpha,digitos))+'))} $' # Errores respecto a rssi(d) xi0 = grssi0(yi) dxi0 = xi - xi0 dxi0mean = np.mean(np.abs(dxi0)) dxi0std = np.std(dxi0, dtype=np.float64) unaecuacion = {'alpha' : alpha, 'beta' : beta, 'eq_latex': fdtxt0, 'intervalox' : [np.min(xi),np.max(xi)], 'error_medio': dyi0mean, 'error_std' : dyi0std, 'eqg_latex' : grtxt0, 'intervaloy' : [np.min(yi),np.max(yi)], 'errorx_medio': dxi0mean, 'errorx_std' : dxi0std, } return(unaecuacion)
Referencias:
Burden R, Faires J, Burden A, Análisis numérico, Décima Edición 8.1 p370.
Chapra C, Canale R. Métodos numéricos para ingenieros, Quinta edición 17.1.2 p469.Numpy org. Least Square function. numpy.linalg.lstsq
Mínimos cuadrados. https://es.wikipedia.org/wiki/M%C3%ADnimos_cuadrados
