5

I'd like to calculate the receiver GPS location myself, without relying on high-level API. A few questions there: I am looking at an Android app that can show GPS satellites, on a diagram like this

enter image description here

However Android (until API 25) does not return satellite pseudorange, which every source I read is required to compute a location yourself. Android API does expose signal-to-noise ratio, elevation, azimuth, and PNR (pseudorandom number) for each satellite in view. Now I believe the elevation and azimuth are given in relation to the receiver (hence the image above), so I was wondering, given this information and precise location of each satellite (the PRN uniquely identifies a satellite AFAIK, I could look that up from an almanac, or ephemerisis file), couldn't we compute a rough location?

This gentleman here can compute the location of a GPS satellite given receivers location, satellite elevation and azimuth. I want to merely do the reverse, given satellite location (from almanac), azimuth and elevation, compute receiver's location.

Another link on azimuth / elevation and Android.

And this person seems to be talking about the distance computation, and even computing the receiver location..?

  • Why not use current receiver location (in developers.google.com/android/reference/com/google/android/gms/…) which will use the GPS receiver as intended? What are you actually trying to do? – BradHards Dec 5 '16 at 10:08
  • I am trying to compute the receiver location myself - I noticed it takes a lot of satellites in view for Android API to compute a lat/lon; when GPS references indicate 4 satellites is enough. Apparently pseudorange is needed for that but Android won't give me raw GPS data until API 25, so I wondered if I could approximate pseudorange from elevation, azimuth, PRN and almanac (ephemerisis) and do the calculation myself. – BBDynSys Dec 5 '16 at 10:24
  • 1
    Since you seem to want to do this as an exercise, the general approach is to create a set of simultaneous linear equations projecting a line from each satellite and then solve the system of equations for a point of intersection using e.g., a least-squares error minimization. Given the precision available you won't be able to do very well but it might be fun to see how well in practice. – J Hayward Dec 9 '16 at 21:50
  • 1
    Calculating the receiver position from satellite elevation and azimuth angles is relatively "simple" geometry and is entirely unrelated to what the GNSS receiver does with pseudoranges etc. The Azimuth and elevation for satellites are calculated by the GNSS receiver based on the satellite positions and the receiver position. The Android app does not calculate a position - this is supplied by the receiver along with the Az and El for the satellites, etc. Do a Google search for "uBlox Neo-6 Protocol Specification" and download that to see what a typical receiver supplies. – Trams Apr 3 '17 at 8:07
  • Further to what @Trams wrote - the receiver in your phone does not measure az or el. It only measures pseudorange. Az and el are calculated after the receiver already has its position. So just use that! – BradHards May 27 '17 at 22:48
1

Here is how to go in the other direction, computing azimuth, elevation based on known receiver and satellite locations (latitude, longitude, altitude). If this can be done, then given many azimuth and elevation data points (and given satellite location) one could solve for unknown receiver location using nonlinear optimization.

The best explanation on how to compute azimuth and elevation is given here. The implementation shown is in Pascal, but the pyorbital package wonderfully converted this code into Python. I pulled out the relevant code and am including it here:

import datetime
import numpy as np

F = 1 / 298.257223563 # Earth flattening WGS-84
MFACTOR = 7.292115E-5 
EPS_COS = 1.5e-12
F = 1 / 298.257223563  # Earth flattening WGS-84
A = 6378.137  # WGS84 Equatorial radius

def jdays2000(utc_time):
    """Get the days since year 2000.
    """
    return _days(utc_time - datetime(2000, 1, 1, 12, 0))


def jdays(utc_time):
    """Get the julian day of *utc_time*.
    """
    return jdays2000(utc_time) + 2451545

def _fdays(dt):
    """Get the days (floating point) from *d_t*.
    """
    return (dt.days +
            (dt.seconds +
             dt.microseconds / (1000000.0)) / (24 * 3600.0))

_vdays = np.vectorize(_fdays)

def _days(dt):
    """Get the days (floating point) from *d_t*.
    """
    try:
        return _fdays(dt)
    except AttributeError:
        return _vdays(dt)

def gmst(utc_time):
    """Greenwich mean sidereal utc_time, in radians.

    As defined in the AIAA 2006 implementation:
    http://www.celestrak.com/publications/AIAA/2006-6753/
    """
    ut1 = jdays2000(utc_time) / 36525.0
    theta = 67310.54841 + ut1 * (876600 * 3600 + 8640184.812866 + ut1 *
                                 (0.093104 - ut1 * 6.2 * 10e-6))
    return np.deg2rad(theta / 240.0) % (2 * np.pi)


def observer_position(time, lon, lat, alt):
    """Calculate observer ECI position.

    http://celestrak.com/columns/v02n03/
    """

    lon = np.deg2rad(lon)
    lat = np.deg2rad(lat)

    theta = (gmst(time) + lon) % (2 * np.pi)
    c = 1 / np.sqrt(1 + F * (F - 2) * np.sin(lat)**2)
    sq = c * (1 - F)**2

    achcp = (A * c + alt) * np.cos(lat)
    x = achcp * np.cos(theta)  # kilometers
    y = achcp * np.sin(theta)
    z = (A * sq + alt) * np.sin(lat)

    vx = -MFACTOR*y  # kilometers/second
    vy = MFACTOR*x
    vz = 0
    return (x, y, z), (vx, vy, vz)   

def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt):
    """Calculate observers look angle to a satellite.
    http://celestrak.com/columns/v02n02/
    utc_time: Observation time (datetime object)
    lon: Longitude of observer position on ground
    lat: Latitude of observer position on ground
    alt: Altitude above sea-level (geoid) of observer position on ground
    Return: (Azimuth, Elevation)
    """
    (pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = observer_position(
        utc_time, sat_lon, sat_lat, sat_alt)

    (opos_x, opos_y, opos_z), (ovel_x, ovel_y, ovel_z) = \
        observer_position(utc_time, lon, lat, alt)

    lon = np.deg2rad(lon)
    lat = np.deg2rad(lat)

    theta = (gmst(utc_time) + lon) % (2 * np.pi)

    rx = pos_x - opos_x
    ry = pos_y - opos_y
    rz = pos_z - opos_z

    sin_lat = np.sin(lat)
    cos_lat = np.cos(lat)
    sin_theta = np.sin(theta)
    cos_theta = np.cos(theta)

    top_s = sin_lat * cos_theta * rx + \
        sin_lat * sin_theta * ry - cos_lat * rz
    top_e = -sin_theta * rx + cos_theta * ry
    top_z = cos_lat * cos_theta * rx + \
        cos_lat * sin_theta * ry + sin_lat * rz

    az_ = np.arctan(-top_e / top_s)

    az_ = np.where(top_s > 0, az_ + np.pi, az_)
    az_ = np.where(az_ < 0, az_ + 2 * np.pi, az_)

    rg_ = np.sqrt(rx * rx + ry * ry + rz * rz)
    el_ = np.arcsin(top_z / rg_)

    return np.rad2deg(az_), np.rad2deg(el_)

On Android using this package, I recorded some satellite values, one example is 22:0.0:331.0:7.0, the values are PRN:SNR:azimuth:elevation. For that PRN, I took the position of the satellite from here (updated regularly). I calculate,

millis = 1492964228447
from datetime import datetime
dt = datetime.fromtimestamp(millis/1000.0)

sat_lon = -109.861429
sat_lat = -6.725577
sat_alt = 20*1000*1000
lon = 13.442383333333332
lat = 52.483086666666665
alt = 0
az,el = get_observer_look(sat_lon, sat_lat, sat_alt, dt, lon, lat, alt)
print az, el  

The result was

293.517066155 -25.1656865911

Not sure what to do with negative elevation but the azimuth seemed close to recorded values on Android. I checked few others, they were okay.

Going in the other direction would be an optimization problem. But of course cleaner solutions are probably possible.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.