Too many gateways leading to too many hyperbolas and LSE ends up at GW[c] itself!!

mkdthanga commented 3 years ago

GW_vs_Tx When a device reached too many gateways (10), then there is too many hyperbolas intersection and as a result, the solved position would be the starting point of the equation which is the GW[c], that is the closest gateway would be solved as Tx position! Is there any solution to avoid this behavior of the LSE resolver? Of course, the Tx is somewhere near the GW[c], but it would be better if it is not exactly showing the gateway location rather somewhere between GW[c] and the next nearest gateway!

Thank you in advance!

PFA screenshot!

Regards, Thangaraj

jurasofish commented 3 years ago

Hard to say why that's happening. looks like you're using code that's different, based on this repo? A couple things you could check:

In your plot it looks like most of the time there are two gateways in the same spot. The labels look like they have two sets of numbers in them. If this is the case you'll probably get some very large or very small numbers creeping in -> numerical instability
One hyperbola looks like a straight line to me. Consistent with previous point of numerical instability.
If you're introducing errors to the arrival timestamps that could cause the problem. Maybe try using fresh random data.

But yeah without the full code and data there's not much I can offer you.

mkdthanga commented 3 years ago

Here is the code for my LoRaWAN Geolocation Calculation:

import numpy as np
import pandas as pd
from scipy.optimize import least_squares
from scipy.optimize import minimize
from pyproj import Proj
from numpy import linalg
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
from multilaterate import get_loci
# UTM coordinates converter
myProj = Proj("+proj=utm +zone=33 +north +ellps=WGS84")

# true location (just a guess, still not sure where the device is! )
tx_gps = [47.8387297, 12.1499166]
tx_utm = myProj(tx_gps[1], tx_gps[0])
# Hyperbola plot parameters
plotwidth= 50000
plot_lines_between_GWs = False
plot_trilateration_circles = False
delta_d = int(100)
tx_square_side = 5e3
rx_square_side = 25e3
rec_time_noise_stdd = 1e-6
max_d = int(20e3)
# Speed of transmission propogation
v = 3e8
# Resolve the tdoa position given the arrival timestamps & gateway locations
def resolve_tdoa(df):
    TOA= np.asarray(pd.to_numeric(df['toa'])).astype(float)
    lat = np.array(pd.to_numeric(df['gateway_latitude'])).astype(float)
    lon = np.array(pd.to_numeric(df['gateway_longitude'])).astype(float)
    num_TOAs = len(TOA)
    GWs = np.zeros((num_TOAs, 2))
    # convert lat, lon to x,y
    for i in range(GWs.shape[0]):
        GWs[i][0], GWs[i][1] = myProj(lon[i], lat[i])
    #guess= GWs.mean(axis=0)
    c = np.argmin(TOA)
   # assign initial guess as nearest gateway
    guess = GWs[c]
    # Call get_loci function
    # loci = get_loci(TOA, GWs, v, delta_d, max_d)
    p_c = np.expand_dims(GWs[c], axis=0)
    t_c = TOA[c]
    # Remove the c GW to allow for vectorization.
    all_p_i = np.delete(GWs, c, axis=0)
    all_t_i = np.delete(TOA, c, axis=0)
    # Position Calculation of Sensor
    def eval_solution(x):
        return (np.linalg.norm(x - p_c, axis=1) - np.linalg.norm(x - all_p_i, axis=1) + v * (all_t_i - t_c))
    res = least_squares(eval_solution, guess)
    return res

# function for hyperbola plots
def plot_hyperbola(df):
    # extract the subset/ individual components of dataframe being grouped
    toa = np.asarray(pd.to_numeric(df['toa'])).astype(float)
    lat = np.array(pd.to_numeric(df['gateway_latitude'])).astype(float)
    lon = np.array(pd.to_numeric(df['gateway_longitude'])).astype(float)
    num_GWs = len(toa)
    GWs = np.zeros((num_GWs, 2))
    TOA = np.zeros(num_GWs)
    # convert lat, lon to x,y
    for i in range(GWs.shape[0]):
        GWs[i][0], GWs[i][1] = myProj(lon[i], lat[i])
        TOA[i] = toa[i]
    # Get the loci.
    loci = get_loci(TOA, GWs, v, delta_d, max_d)
    # Plot GWs and transmission location.
    fig, ax = plt.subplots(figsize=(5, 5))
    plt.title('Given GWs & Tx Location')
    plt.xlabel('UTMx (longitude in meters)')
    plt.ylabel('UTMy (latitude in meters)')
    #print('num of geoloc frames:', GWs.shape[0])
    for i in range(GWs.shape[0]):
        x = GWs[i][0]
        y = GWs[i][1]
        ax.scatter(x, y)
        ax.annotate('GW ' + str(i), (x, y))
    ax.scatter(tx_utm[0], tx_utm[1])
    ax.annotate('Tx', (tx_utm[0], tx_utm[1]))
    for locus in loci:
        # print('loci:' + str(loci)+ ':', locus[0], locus[1])
        ax.plot(locus[0], locus[1])
    plt.show()
    # plot the same hyperbola using the derived expressions.
    fig, ax = plt.subplots(figsize=(5, 5))
    plt.title('Hyperbola using the derived expressions')
    plt.xlabel('UTMx (longitude in meters)')
    plt.ylabel('UTMy (latitude in meters)')
    c = np.argmin(TOA)  # GW that received first.
    p_c = GWs[c]
    t_c = TOA[c]
    x = np.linspace(GWs[i][0] - plotwidth, GWs[i][0] + plotwidth, 100)
    y = np.linspace(GWs[i][1] - plotwidth, GWs[i][1] + plotwidth, 100)
    x, y = np.meshgrid(x, y)
    for i in range(num_GWs):
        if i == c:
            continue
        p_i = GWs[i]
        t_i = TOA[i]
        plt.contour(x, y,(np.sqrt((x - p_c[0]) ** 2 + (y - p_c[1]) ** 2)-np.sqrt((x - p_i[0]) ** 2 + (y - p_i[1]) ** 2)+ v * (t_i - t_c)), [0])
    # Solve the location of the Sensor.
    c = np.argmin(TOA)
    p_c = np.expand_dims(GWs[c], axis=0)
    guess= GWs[c]
    t_c = TOA[c]
    # Remove the c GW to allow for vectorization.
    all_p_i = np.delete(GWs, c, axis=0)
    all_t_i = np.delete(TOA, c, axis=0)
    # Position Calculation of Sensor
    def eval_solution(x):
        """ x is 2 element array of x, y of the Sensor"""
        return (np.linalg.norm(x-p_c,axis=1)-np.linalg.norm(x-all_p_i,axis=1)+v*(all_t_i-t_c))
    # Find a value of x such that eval_solution is minimized.
    # Remember the receive times have error added to them: rec_time_noise_stdd.
    res = least_squares(eval_solution, guess)
    # And now plot the solution.
    fig, ax = plt.subplots(figsize=(5, 5))
    plt.title('Solution using LSE')
    plt.xlabel('UTMx (longitude in meters)')
    plt.ylabel('UTMy (latitude in meters)')
    for locus in loci:
        ax.plot(locus[0], locus[1])
    for locus in loci:
        ax.plot(locus[0], locus[1])
    ax.scatter(tx_utm[0], tx_utm[1], color='red')
    ax.annotate('Actual', (tx_utm[0], tx_utm[1]))
    ax.scatter(res.x[0], res.x[1], color='blue')
    ax.annotate('Solved', (res.x[0], res.x[1]))
    plt.show()
    return res

# extract the geoloc packets from .csv file!
lora_geoloc_packets = 'lora_geoloc_packets.csv'

df = pd.read_csv(lora_geoloc_packets, sep=',', encoding='latin-1', header= 'infer', engine='python')

print('df:', df)

# concatenate toa seconds and toa nanoseconds
df['toa'] = df['toa_sec'].astype(str).str.cat(df['toa_nsec'].astype(str), sep='.').copy()
grouped = df.groupby(df.frame_cnt)
df = grouped.get_group(df.frame_cnt.unique()[-1])
print('df:', df)

# call the above functions
resolve_tdoa(df)      # returns the numerical tdoa result only 
plot_hyperbola(df)   # returns the hyperbola plots

jurasofish / multilateration

Too many gateways leading to too many hyperbolas and LSE ends up at GW[c] itself!! #3