Convert Latitude, Longitude, Altitude to North, East, Down and Find Distance Between Base and Rover

In order to achieve accurate and precise positioning on a moving platform like a drone, an RTK (real time kinematics) GNSS is the best option. The u-blox C94M8P achieves RTK positioning by using two GNSS receivers: the rover mounted to the drone and the base stationary on land.

For many of the imaging tasks we wanted to perform in my undergrad research group, we needed to convert the typical latitude, longitude, altitude (LLA) coordinate system to a north, east, down (NED) coordinate system.

For more information about what NED is and why it’s often used in aerospace, this Wikipedia article is a good start.

The short version is NED measures the north, east, and down directions from an origin of your choosing. It’s a convenient system to use if you want to know where a moving platform is relative to some reference point, as is the case with our drone’s rover-base positioning system. The drone’s initial position can also be used as the origin, but in our application with two GPS receivers, the base is a natural choice.

ned

An advantage of NED is you can easily obtain distance between base and rover by simply implementing the 3D distance formula:

\mathbf{distance=\sqrt{north^{2}+east^{2}+down^{2}}}

While MATLAB has a built-in function, geodetic2ned, which converts geodetic (LLA) coordinates to NED, it’s only available in their mapping toolbox, which requires an additional purchase. So, I did some research and created my own function, which I named LLAtoNED. I share it freely in the name of open source.


function [time, north, east, down, dist] = LLAtoNED(baseFile, roverFile)

%{
    File:           LLAtoNED.m
    Author:         Brandi Downs
    Last Updated:   12.11.2018

    LLAtoNED.m reads in a file specified by the user, looks for the
    parameters given below, and stores them to a cell array for further
    processing.

    This script currently stores data only from $GNGGA sentences.

    This script assumes the base GPS unit is stationary.

    Parameters:
    TIME
    LAT
    LONG
    FIX
    HDOP
    ALT
    SEP
%}

%%  Open base and rover txt files, read in NMEA sentence data, store data in new arrays
roverData = fopen(roverFile,'r');      % open rover file for reading
row = 1;        % begin with row 1

% initialize a matrix for each parameter you want to store
gngga = []; % stores all data from GNGGA sentences
time =  []; % time in UTC
lat =   []; % latitude
long =  []; % longitude
fix =   []; % fix type or quality indicator (see note above)
sats =  []; % satellites used
hdop =  []; % horizontal dilution of precision
alt =   []; % altitude
sep =   []; % geoidal separation

format long;

while ~feof(roverData)    % while not end of file

    line = fgetl(roverData);  % get next line
    line = regexp(line,',','split');    % split line into cell array according to commas

    if ~any(strcmpi('$GNGGA',line))     % if line doesn't contain 'GNGGA', continue to next iteration of loop
        continue
    end

    for i = 2:length(line)
        gngga(row,i-1) = str2double(line{i});   % store all data to gngga matrix

        % store time parameter to column vector
        time(row,1) = str2double(line{2});

        % convert degree, minute, decimal minute to decimal degrees
        lat(row,1) = str2double(line{3})/100;   % move decimal place left 2
        latdeg = floor(lat(row,1));             % extract integer part
        latdec = lat(row,1) - latdeg;           % extract decimal part
        latdec = latdec/60;                     % convert decimal minutes to decimal degrees
        latdec = latdec*100;                    % move decimal place right 2
        lat(row,1) = latdeg + latdec;           % add integer and decimal to complete calculation
        % do same for longitude
        long(row,1) = str2double(line{5})/100;  % move decimal place legt 2
        longdeg = floor(long(row,1));           % extract integer part
        longdec = long(row,1) - longdeg;        % extract decimal part
        longdec = longdec/60;                   % convert decimal minutes to decimal degrees
        longdec = longdec*100;                  % move decimal place right 2
        long(row,1) = longdeg + longdec;        % add integer and decimal to complete calculation
        % continue to store remaining parameters
        fix(row,1) = str2double(line{7});
        sats(row,1) = str2double(line{8});
        hdop(row,1) = str2double(line{9});
        alt(row,1) = str2double(line{10});
        sep(row,1) = str2double(line{12});
    end

    row = row + 1;  % next row
end

fclose(roverData);

% Do same for base data
% Only keep lat, long, alt
baseData = fopen(baseFile,'r');     % open base file for reading
row = 1;
latBase =  []; % latitude
longBase = []; % longitude
altBase =  []; % altitude
while ~feof(baseData)   % while not end of file
    line = fgetl(baseData); % get next line
    line = regexp(line,',','split');    % split line into cell array according to commas
    if ~any(strcmpi('$GNGGA',line))     % if line doesn't contain 'GNGGA', continue to next iteration of loop
        continue
    end

    for i = 2:length(line)
        % get base latitude and convert to decimal degrees
        latBase(row,1) = str2double(line{3})/100;
        latBasedeg = floor(latBase(row,1));
        latBasedec = latBase(row,1) - latBasedeg;
        latBasedec = latBasedec/60;
        latBasedec = latBasedec*100;
        latBase(row,1) = latBasedeg + latBasedec;
        % get base longitude and convert to decimal degrees
        longBase(row,1) = str2double(line{5})/100;
        longBasedeg = floor(longBase(row,1));
        longBasedec = longBase(row,1) - longBasedeg;
        longBasedec = longBasedec/60;
        longBasedec = longBasedec*100;
        longBase(row,1) = longBasedeg + longBasedec;
        % get base altitude
        altBase(row,1) = str2double(line{10});
    end
    row = row + 1;  % next row
end
fclose(baseData);

%%  Convert from lat, long, altitude (LLA) to north, east, down (NED) coordinate system

% lat0, long0, alt0 will form the origin of the NED coordinate system
lat0 = median(latBase);
long0 = median(longBase);
alt0 = median(altBase);

% geodetic frame parameters
Rea = 6378137;          % semi-major axis (m)
f = 1/298.257223563;    % flattening factor
Reb = Rea*(1-f);        % semi-minor axis
ecc = sqrt(Rea^2 - Reb^2)/Rea;      % first eccentricity
Ne = zeros(1,numel(lat));      % prime vertical radius of curvature
for k = 1:numel(lat)
    Ne(k) = Rea./sqrt(1-(ecc^2)*(sind(lat(k))).^2);
end
Neref = Rea./sqrt(1-(ecc^2)*(sind(lat0)).^2);   % prime vertical radius of curvature for reference point

%% Geodetic to Earth-Centered Earth-Fixed (ECEF) -- Intermediate step

% preallocate arrays for speed
xe = zeros(1,numel(lat));
ye = zeros(1,numel(lat));
ze = zeros(1,numel(lat));

% convert rover lat, long, alt (LLA) to ECEF
for k = 1:numel(lat)
    xe(k) = (Ne(k) + alt(k))*cosd(lat(k))*cosd(long(k));
    ye(k) = (Ne(k) + alt(k))*cosd(lat(k))*sind(long(k));
    ze(k) = (Ne(k)*(1-ecc^2)+alt(k))*sind(lat(k));
end

% convert reference position (origin in NED frame) from LLA to ECEF
xeref = (Neref + alt0)*cosd(lat0)*cosd(long0);
yeref = (Neref + alt0)*cosd(lat0)*sind(long0);
zeref = (Neref*(1-ecc^2)+alt0)*sind(lat0);

Pe = [xe;ye;ze];
Peref = [xeref;yeref;zeref];    

%% ECEF to NED

% Rne is the rotation matrix from ECEF frame to local NED frame
Rne = [-sind(lat0)*cosd(long0)  -sind(lat0)*sind(long0) cosd(lat0); ...
       -sind(long0)             cosd(long0)             0; ...
       -cosd(lat0)*cosd(long0)  -cosd(lat0)*sind(long0) -sind(lat0)];

% convert ECEF to local NED frame
Pn = zeros(5,numel(lat));
for k = 1:numel(lat)
    Pn(1,k) = time(k);
    Pn(2:4,k) = Rne*(Pe(:,k) - Peref);
    Pn(5,k) = sqrt((Pn(2,k))^2 + (Pn(3,k))^2 + (Pn(4,k))^2);
end

% use num2str(x,'%4.10f') to display full numbers
time = Pn(1,:);
north = Pn(2,:);
east = Pn(3,:);
down = Pn(4,:);
dist = Pn(5,:);

LLAtoNED.m is a MATLAB function and hence must be called as a function command. By typing:

[time,north,east,down,dist] = LLAtoNED(‘nameOfBaseFile.txt’,’nameOfRoverFile.txt’)

in the command window, LLAtoNED will read in text files for the base and rover position data, do the coordinate conversions from LLA to NED, and output arrays for time, north, east, down, and distance between base and rover at each observation data point. For the ublox C94M8P in wireless RTK mode, the observation rate — or how many times the base and rover record their position data — is a maximum of 3 Hz max (3 times per second).

Enjoy and please leave comments if you have any questions or feedback!

Note: For the coordinate conversion matrices, I used this very helpful reference document: Coordinate Systems and Transformations

How gpsd can interfere with your GNSS receiver and how to fix it

Pi and u-blox

Raspberry Pi 3 Model B and u-blox C94-M8P

In order to work with the u-blox C94-M8P on the Raspberry Pi, I installed a GPS interface called gpsd on the Pi. Documentation on gpsd can be found here. And this is the tutorial I used to install it on the Pi. However, there is one issue with gpsd that neither tutorial mentioned.

It turns out, when gpsd attempts to reconfigure certain Bluetooth- or USB-connected GPS receivers, it can interfere in strange ways with their functioning. I found that certain gpsd commands interfered with the C94-M8P’s ability to reach RTK mode. The receiver would initially reach RTK mode, then quickly lose it and achieve only a 2D/3D fix while gpsd was running.

After much digging, I found a solution. You can configure gpsd to run with the -b option, which is a read-only mode that prevents it from writing to the GPS receiver. This can be achieved when you install it, by passing -b to the Options to gpsd prompt, or, if you have already installed gpsd, you simply need to change the default parameters, as detailed below:

  • Start up your Pi and open the default file for gpsd with:

sudo nano /etc/default/gpsd

  • Set the parameters to look like this:

gpsd default settings

  • That’s it. Exit and reboot.

From the gpsd man page:

-b

Broken-device-safety mode, otherwise known as read-only mode. A few bluetooth and USB receivers lock up or become totally inaccessible when probed or reconfigured; see the hardware compatibility list on the GPSD project website for details. This switch prevents gpsd from writing to a receiver. This means that gpsd cannot configure the receiver for optimal performance, but it also means that gpsd cannot break the receiver. A better solution would be for Bluetooth to not be so fragile. A platform independent method to identify serial-over-Bluetooth devices would also be nice.

The read-only mode was a simple solution to a problem that plagued our group for over a week. If you are finding your Pi-connected USB or Bluetooth GPS mysteriously goes in and out of RTK mode and you have exhausted other causes, such as poor satellite reception at the base receiver or incorrect RTCM messages at the roving receiver, try checking your gpsd default parameters and adding this option.

Decoding NMEA Sentences

NMEA sentences (pronounced nee-ma), are a standard format of data output for all GPS receivers. In other words, it’s the language GPS receivers use to communicate the data they produce and receive, such as time, latitude, longitude, altitude, GPS health, speed, etc.

NMEA stands for National Marine Electronics Association, the agency responsible for standardizing the language. The sentences were created as a means for marine electronics to all speak the same language, thus enabling digital communication between different devices.

The sentence structure consists of a string of comma separated values, beginning with a ‘$’ and ending with a checksum. They look like this:

nmea sentences

Each sentence is referred to as a message and tells useful information about the receiver and its positioning. The output above was obtained from a u-blox C94-M8P. Let’s take a closer look at the GNGGA sentence to understand what each field means. Note the contents of each field will be the same for all GGA sentences but minor variations may occur from one GNSS receiver to the next.

NMEA sentence

Message Identifier:  Identifies the source and type of information. The first two letters denote the Talker ID, identifying the source of the information. This sentence is a GN message, or Global Navigation Satellite System (GNSS), meaning a combination of navigation systems was used to obtain the message data. Common talker ID’s are listed below.

Talker ID Constellation Country
GA Galileo European Union
GB BeiDou China
GL GLONASS Russia
GP GPS United States
GN combination (GPS+GLONASS) multiple

The next three letters define the message content. GGA denotes GPS fix data. Some common messages are listed below.

Message ID Meaning
GGA GPS position, time, and fix data
GLL Latitude and longitude data
GNS GNSS position, time, and fix data
GRS GNSS range residuals (error-related data)
GSA Satellite information and dilution of precision (DOP: see HDOP below)

Time:  Universal Time Coordinated (UTC) in hhmmss.ss.

Latitude:  Latitude in decimal degrees (ddmm.mmmm).

North/South Indicator:  North or south of the Equator.

Longitude:  Longitude in decimal degrees (ddmm.mmmm).

East/West Indicator:  East or west of the Prime Meridian.

Quality Indicator:  The quality of the GNSS data — what kind of fix the receiver has obtained.

0 — No fix
1 — Standard GPS (2D/3D) fix
2 — Differential GPS (DGPS) fix
3 — Precise Positioning System (PPS) fix — for government use only
4 — RTK fixed solution
5 — RTK float solution
6 — Estimated (Dead Reckoning, or DR) fix

Satellites Used:  Number of satellites used in solution.

HDOP:  Horizontal Dilution of Precision — a measure of confidence in the solution. Lower numbers indicate a higher confidence level. Thus, the 99.99 shown here denotes highly inaccurate measurements. Refer to Wikipedia’s Dilution of Precision article for more information.

Altitude:  Altitude above mean sea level.

Units:  Units used for altitude or geoidal separation (M = meters).

Geoidal Separation:  Difference between the reference ellipsoid and mean sea level based on the geodetic model WGS84.

DGPS Station ID:  ID of the differential reference station used for DGPS. Blank when DGPS is not used.

Checksum:  Used for data validation.

 

The GNGGA sentence provides all the data one needs for most GNSS analysis and is what we’ll use when we write our script converting latitude, longitude, altitude (LLA) to a north, east, down (NED) coordinate system. However, other sentences may be useful, depending on your navigation needs. Below are some great resources I used in compiling this information.

u-blox 8 / M8 Receiver Description (pdf) — Includes a full breakdown of all NMEA sentences used in u-blox 8 / M8 receivers.

Trimble’s NMEA-0183 Overview

ESRI’s article on Mean Sea Level, GPS, and the Geoid

Eric Raymond’s extensive page on NEMA sentences

u-blox tutorials

For my research project, I’ve worked pretty extensively with a certain Global Navigation Satellite System (GNSS) application board: the u-blox C94-M8P, which includes two u-blox NEO-M8P-2 receivers. It’s a highly accurate GNSS module with tons of cool functionalities (like anti-jamming capabilities) but takes quite the learning curve to be able to understand its various uses and settings.

So I’m starting this website by writing up a series of tutorials explaining some of the things I’ve learned about the C94-M8P and GNSS, such as:

  • how GPS works
  • all about NMEA sentences
  • what are RTCM messages
  • how RTK GPS works
  • using u-center, u-blox’s evaluation software
  • saving configuration settings with u-center
  • configuring the C94-M8P for base/rover mode
  • enabling RTCM messages on the C94-M8P
  • achieving RTK mode on the C94-M8P
  • troubleshooting why the receiver won’t achieve RTK fixed mode
  • using gpsd to interface with the C94-M8P on the Raspberry Pi
  • how gpsd can interfere with a GNSS receiver and how to fix it
  • converting latitude, longitude, and altitude to a north, east, down coordinate system
  • finding distance between base and rover

I’ve gathered this knowledge through a combination of research, reading u-blox manuals, browsing their Q&A forum (and especially learning from clive1), and my own trial and error. My hope is to contribute to the free library of online knowledge that I have benefited from and help others with similar technical problems to those I encountered.