## 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.

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 ## 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:

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.

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
• 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.