gpsd/www/ppp-howto.adoc

= Precise Point Positioning (PPP) HOWTO
Gary E. Miller <gem@rellim.com>
13 Jan 2021
:author: Gary E. Miller
:date: 13 January 2021
:description: This document is a guide getting high accuracy from your GPS using Precise Point Positioning (PPP).
:email: <gem@rellim.com>
:keywords: Precise Point Positioning, PPP, GPSD, GPS
:robots: index,follow
:sectlinks:
:toc: macro

include::inc-menu.adoc[]

== WARNING

This document assumes you are using gpsd versions 3.22 or git
head.  Using older versions will fail in strange ways.

== Introduction

This document is a guide getting high accuracy static positions from
your GPS using Precise Point Positioning (PPP). The rare few that have a
GPS that output raw measurement data for L1 and L2 can achieve absolute
accuracy of around 3 cm.  The u-blox ZED-F9P and ZED-F9T can be better
than 1 cm. The lucky owners of an L1 GPS that outputs raw measurements
can get about 0.5 m. The majority will only be able to get somewhat
better than 1.5 m using simple averaging.

Patience is required.  For best results 6 to 24 hours of data is
required.  Post processing time may double that.

This document is not about getting high precision dynamic positions from
your GPS.  This is about after the fact Post Processed Positions (PPP).
High precision real time dynamic positions requires Real Time Kinematic
(RTK) which will not be discussed here.  RTK users may still want to
read this document.  The RTK Best practice is to determine the position
of the base using PPP before sending the RTK data from the base to its
rovers.

== Background

To compute a location fix a GNSS receiver takes measurements (raw data)
of the relative time of arrival of signals broadcast from the GNSS
satellites.  The receiver also knows the predicted orbit (ephemeris) of
each satellite.  The collection of ephemeris from all the satellites is
known as the ephemerides.  With these two sets of data, the CPU in the
receiver can compute a fix in position in time and space using iterative
techniques.

A major problem is that the predicted orbits are often off by one meter
or more.  Ground stations bounce lasers off the individual satellites as
they pass overhead and use this new data to compute the actual orbits
of the satellites.  Using this new ephemeris data, when it becomes
available, combined with the receiver's raw data, better fixes can be
computed.  This is the basis of PPP.

There are other smaller error sources that PPP can reduce.  The effects of
ocean tides on the Earth's crust,  Clock drift of the satellite clocks.
Troposphere and Ionosphere effects, and more.

The first take on the actual epherides may take 90 minutes to be published.
The ephemerides may be further refined for up to three weeks.

There are two flavors of PPP: static, and kinematic.  Kinematic PPP
simply takes an existing time series of measurements and computes a new
time series of fixes using the actual orbits, instead of the predicted
orbits.  Static PPP assumes all the raw measurements were taken at the
same location over a long duration.  Long is typically 15 minutes to
48 hours.  All the computed fixes are then averaged to result in one
precise location for the receiver's antenna.

Now in 2020, most PPP services only use ephermides from GPS satellites,
with more and more adding GLONASS.  As time passes, more constellations
will be used, leading to slightly increased accuracy.

The earth revolves around its axis every 24 hours, give or take a few
milli seconds.  GPS satellites revolve around the Earth's axis about
every 12 hours, give or take.  PPP will provide best results if your
data is 12, 24 or 36 hours.  This allows some orbital errors to average
out.

== Requirements

This document assumes that you have installed gpsd version 3.22 or later,
basic knowledge of how to use _gpsd_.  Before continuing you should know how
to start and stop _gpsd_, and how to use _cgps_ to see you current
position and fix status.

This document also assumes that you have a GNSS receiver connected to a
local GPSD demon and that _cgps_ shows a stable 3D fix. This will allow
Simple Averaging. You will also need Python and _gnuplot_ installed. The
Python and _gnuplot_ do not need to be installed on the host that is
connected to the GPS, they are merely needed for post processing.

For basic PPP (0.5 m) you will need a GPS that outputs L1 raw
measurement data that _gpsd_ understands. Currently that limits you to a
Javad GPS that support the GREIS language, or a u-blox GPS that support
the UBX-RXM-RAWX messages.

u-blox GPS known to support UBX-RXM-RAWX are: -M8T,-M8F, ZED-F9P and
ZED-F9T.

For advanced PPP (3 cm) you will need a GPS that outputs L1 and L2 raw
measurement data that _gpsd_ understands. Currently that limits you to a
Javad GPS that support the GREIS language or the u-blox 9 series.

Finally, patience is required.

== Results

The end goal of this process is to determine the latitude, longitude and
altitude of your GPS antenna as precisely as possible.  Additionally
the Circular Error of Probability (CEP) will be determined.

CEP, also known as the CEP(50), is the radius of a circle, centered on
the mean, whose boundary includes 50% of the GPS positions.  You probably
do not want to base any navigation or surveying on a 50% probability.

More interesting is the CEP(95) which includes 95% of the measurements.

== ITRF00, WGS84, NAD83 and Ellipsoids

Ever noticed how two "accurate" GPS placed side by side can give wildly
different latitude, longitude, and especially altitude for the same
spot?

All GNSS systems compute positions using ECEF (earth-centered,
earth-fixed) coordinates. After an ECEF position is calculated, it is
converted into latitude and longitude using a datum.  So many to
choose from.  You've probably heard of WGS84, NAD83, and maybe ITRF.

NAD83 is pinned to the North American tectonic plate. WGS84 is pinned
to the Reference Meridian (near the Greenwich Meridian). NAD83 is the
official datum in the USA and Canada, and is used by the FAA.  WGS84 is
the official datum of GPS and the US Department of Defense.

In 1987 the difference between NAD83 and WGS84 was not measurable. Since
then the tectonic plates have moved. In 2018 the two datums can differ
by more than 2 meters in the continental USA.

It is common when using NAD83 to also specify the year (epoch) of the
measurements.  This allows archival, and current, data to be used
to similar accuracy.

It gets worse. Most PPP services support the International Terrestrial
Reference Frame (ITRF). ITRF is pinned to a Celestial Reference Plane
(CRF).

It gets worse. There is not one WGS84, but many: WGS 1984 (ORIG),
WGS84(G730), WGS84(G873), WGS84(G1150), WKID: 4326, and more. There is
not one ITRF but many: ITRF91/92, ITRF94/96, ITRF00, ITRF08 and more.

Original ITRF and WGS84 differ by less than 1 meter, which is huge
for the purposes here. ITRF2014 and WGS84(G1762) differ by a few
centimeters.

It gets worse.  Two expensive GPS often differ in altitude by over 60
feet.  The Earth is not a perfect sphere. It is more pear shaped.  GPS
approximate this with an ellipsoid, usually some version WGS84.  Then
altitude is calculated as height above the ellipsoid (HAE).

Most people do not consider the altitude as the height above the
ellipsoid, but as the height above Mean Sea Level.  MSL is the same
as pressure altitude (when corrected for temperature and barometric
pressure), but different from HAE.

Mean Sea Level has had nothing to with the level of the sea, it relates
to some stakes driven into the ground 100 hundred years ago that seemed
at the time to be roughly mean sea level at that point.

Pressure altitude is not some sort of absolute geometric altitude,
it is related to the gravity under the position being measured.  So
local gravity affects local MSL, but not local HAE.  If the measured
position is over a high density rock, like iron ore, then the gravity is
higher, and the pressure is increased over the simple HAE.  Conversely
ocean water is less dense and has the opposite effect.  This is very
noticeable in the Hawaiian Islands.

Fun fact: gravity is at its maximum at 'sea level'.

Many different datums can be used to calculate height. These datums
are based on different ellipsoids used to approximate sea level. The
two used by CSRS-PPP are CGVD2013 and CGDV28(HT2_0). GCVD2013 is the
standard datum in Canada since 2013.  The standard in the USA is the
North American Vertical Datum of 1988 (NAVD 88). Many GPS use the WGS84
Ellipsoid as the vertical datum. The WGS84 Ellipsoid is from the same
organization as the WGS84 coordinate system, but not part of WGS84.

More refined datums, like the World Gravity Model, WGM2012, also take
into account more gravitation effects.

Clearly there is no point knowing your precise position to a few cm
if you are not certain of your datum and vertical datum (ellipsoid),
with epoch. This will be important later as you are asked to input your
choice of horizontal and vertical datums to your PPP service.

== Averaging

The first technique covered, Simple Averaging, works with any GPS that
is supported by _gpsd_.  For best results a minimum of 6 hours, and
preferably 24 hours, of continuous observations are required.

_gpsprof_ will be used to gather 24 hours of position data and then
output a plot file. The plot file is fed into _gnuplot_ to turn it
into a png image file. The image will contain a scatter plot of all
the positions reported by your GPS, as well as summary statistics. The
statistics include the mean latitude, mean longitude, mean altitude and
other computed values.

The procedure is simple:

. Verify your GPS is communicating with _gpsd_ by running _cgps_ and
confirming that you have a stable 3D fix.

. Collect 24 hours of data in a plot file: `gpsprof -n 2880 -T pngcairo > scatter.plot`

. Convert the plot to a png: `gnuplot < scatter.plot > scatter.png`

. Display the png with your favorite image viewer.  To use _display_
from _Imagemagick_: `display scatter.png`

There are many possible adjustments to the above procedure.

Maybe you want to collect just 10 minutes of data (20 epochs at 30 second
interval) to verify that your
tool-chain is working before doing a 24 hour run. Simple, just change
`gpsprof -n 2880` to `gpsprof -n 20` and then proceed as above.

Maybe your _gpsd_ host does not have Python installed.  Just run _gpsprof_
remotely.  On the host you will need to run _gpsd_ with the `-g` parameter so
that it can be accessed over the network.  Then run _gpsprof_ on a
remote host that supports Python this way:
`gpsprof -n 2880 -T pngcairo [hostname] > scatter.plot`

Depending on your GPS, your GPS antenna, and your sky view, you may get
a CEP(95) of around 1.5 m.

== Precise Point Positioning (PPP)

Plain GPS determine their position by measuring the distance to several
GPS satellites and calculating a position solution. The main limitation
is that the position of any GPS satellite is not known to better than a
meter or two in real time.

PPP uses the raw GPS measurements from a worldwide network of precisely
fixed ground receivers to precisely calculate the actual orbits of
all the satellites. "Ultra Rapid" orbits take about 90 minutes to be
available. "Rapid" orbits take a day. The most accurate orbits ("Final")
take around 14 days to determine.

To use these orbits you need to collect the raw measurements from your
GPS, then upload them to a service to compute a more precise fix.
Receiver Independent Exchange Format (RINEX) files are the standard
for sending your raw measurement data.  _gpsd_ uses RINEX Version 3
(RINEX 3).

Most PPP services have many limitations making them unsuitable for
our purposes.  Some limitations include: open only to paid subscribers,
require L1 and L2 raw data, and/or use proprietary data formats.

There is one online service that is free to all (requires registration),
accepts L1 only raw data, and accepts RINEX 3 files: Natural Resources
Canada (NRCAN).  Their tool is at https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php

Trimble has a free to all (requires registration) service that requires
L1 and L2 observations in RINEX 3.
  Their
tool is at: https://trimblertx.com/Home.aspx

== PPP Configuration

Before you can collect raw data from you GPS, you must configure it to
output raw data.  This configuration will not be the default configuration
that _gpsd_ applies to your GPS by default.

The raw data can be quite large, so be sure your GPS serial port speed
is set to 57,600, or higher.

Many of the configuration steps are order dependent. If in doubt, start
over from the beginning. Be sure that _gpsd_ is running and that _cgps_
shows that you have a stable 3D fix.

=== u-blox configuration

This section is only for u-blox users.

Be sure your serial port speed is high enough:

...................................
$ gpsctl -s 115200
...................................

Disable all NMEA messages, and enable binary messages:

...................................
$ ubxtool -d NMEA
$ ubxtool -e BINARY
...................................

To start simple, disable all constellations, except GPS (and QZSS):

...................................
$ ubxtool -d GLONASS
$ ubxtool -d BEIDOU
$ ubxtool -d GALILEO
$ ubxtool -d SBAS
$ ubxtool -e GPS
...................................

Verify that only GPS and QZSS are enabled.  Otherwise the u-blox 8 will
not output raw measurement data.  You may enable the other constellations
with a u-blox 9, but support for non-GPS in PPP services is limited.

...................................
$ ubxtool -p CFG-GNSS
[...]
UBX-CFG-GNSS:
 Ver: 0 ChHw; 20 ChUse: 20, Blocks: 7
 gnssId: GPS TrkCh: 8 maxTrCh: 16, Flags: 0x1 01 00 01
  L1C/A enabled
 gnssId: SBAS TrkCh: 1 maxTrCh: 3, Flags: 0x1 01 00 00
  L1C/A
 gnssId: Galileo TrkCh: 4 maxTrCh: 8, Flags: 0x1 01 00 00
  E1OS
 gnssId: BeiDou TrkCh: 8 maxTrCh: 16, Flags: 0x1 01 00 00
  B1I
 gnssId: IMES TrkCh: 0 maxTrCh: 8, Flags: 0x3 01 00 00
  L1
 gnssId: QZSS TrkCh: 0 maxTrCh: 3, Flags: 0x5 01 00 01
  L1C/A enabled
 gnssId: GLONASS TrkCh: 8 maxTrCh: 14, Flags: 0x1 01 00 00
  L1OF
[...]
...................................

Enable the good stuff, the raw measurement messages:

...................................
$ ubxtool -e RAWX
...................................

Verify raw data messages are being sent:

...................................
$ ubxtool | fgrep RAWX
...................................

You should see this output that confirms you are seeing raw measurement
data from the GPS:

...................................
UBX-RXM-RAWX:
UBX-RXM-RAWX:
...................................

After you have completed these steps, do not restart _gpsd_.  If you restart
_gpsd_ then you must restart the configuration from the beginning.

=== Javad (GREIS) configuration

The section is only for users of Javad GPS supporting the GREIS
language.

Be sure your serial port speed is high enough.  use _zerk_, _gpsctl_
may be flaky:

...................................
$ zerk -S 115200
...................................

Disable all messages, then enable raw data messages:

...................................
$ zerk -p DM
$ zerk -e RAW
...................................

GREIS will happily send data for all satellites seen, but PPP services
only use GPS and maybe GLONASS. Disable all constellations, except GPS
and QZSS:

...................................
$ zerk -d COMPASS
$ zerk -d GALILEO
$ zerk -d SBAS
$ zerk -e GPS
...................................

Verify that only GPS and QZSS are enabled:

...................................
$ zerk -p CONS
zerk: poll CONS
RE: %cons%/par/pos/sys={gps=y,glo=y,gal=n,sbas=n,qzss=n,comp=n,irnss=n}
...................................

Verify raw data messages are being sent:

...................................
$ zerk -v 2 | fgrep '[PC]'
...................................

You should see this output that confirms you are seeing raw measurement
data from the GPS:

...................................
[PC] cp 199266957.2307 113917941.9777 122453730.9966 108761050.8140 105892190.3611 199725013.5654 117456220.7611 125484683.4227 199977132.8627 126963987.0936 121945102.6244 114688862.4874 140928054.2405 128350477.4361 129924383.6416 199424925.2522 126077127.2204 126780423.4782 120799412.3999
[PC] cp 199266051.1359 113915242.3018 122452018.0540 108761104.8641 105890706.6420 199724109.4819 117454519.9705 125481341.1019 199976227.8647 126966862.6124 121942821.9832 114690162.3442 140924407.3081 128351475.5908 129920370.5866 199424017.5063 126073289.2387 126782833.2288 120800324.7775
...................................

After you have completed these steps, do not retart _gpsd_.  If you restart
_gpsd_ then you must restart the configuration from the beginning.

== Acquire the Raw Data

Configuration complete. Collect 24 hours of samples at 30 second
intervals, save the raw data as RINEX 3 format in the file _today.obs_.
Collecting data at a rate faster than 30 second intervals may degrade
your results.  Trimble will average data to 10 second intervals if
the data rate is faster than 10 seconds. Start the long process:

...................................
$ gpsrinex -i 30 -n 2880 -f today.obs
...................................

Now is a good time to go the NRCAN's CSRS-PPP page and sign up
for a free account.  You need this account to be able to upload the
RINEX 3 file _today.obs_ to their free PPP service for processing.
https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php

Take a break. You now have 24 hours to contemplate the answer to the
ultimate question of life, the universe, and everything.

== Post Process the Raw Data

More waiting.  Before you can post process your data, the PPP service
must be ready for it.  Depending on the service it can take from 10 to
60 minutes before you can upload your new data.  For best results you
should wait 2 weeks.

The following two services are known to work with _gpsrinex_. CSRS-PPP
will accept L1 only data, trimble RTX requires L1 and L2 data.  Try
both, with the same data set, if you can.  That will show you that their
sigma's are "optimistic".

=== CSRS-PPP

After _gpsrinex_ is complete, you need to login to CSRS-PPP and
upload the RINEX 3 file.  After login you will be taken to the upload
page.  Enter your email address, so the results can be emailed to you.

Select processing mode of Static, using the ITRF datum.  Use the "Browse"
button to select the _today.obs_ file with your raw observations.  Then
push "Submit to PPP".

All done, except for more waiting.  You will receive an email from NRCAN
maybe within minutes, maybe up to 36 hours later, with a link to a file
called: full_output.zip.  Unzip, and Voila!  Inside is a PDF file with
your precise position, and other goodies.

=== Trimble RTX

Before uploading today.obs to Trimble you will need to change the _.obs_
extension to _.YYo_, where YY is the 2-digit year.  Then proceed as
above with CSRS-PPP.

=== GAPS

The University of New Brunswich has an online PPP service.  They call
it GNSS Analysis and Positioning Software (GAPS).  GAPS requires
observations from the L2 P signal or L5 signal.  No u-blox chip
follows the L2 P signal.

== References

Wikipedia has a little information on PPP:
https://en.wikipedia.org/wiki/Precise_Point_Positioning

Information on how different datums differ:
https://confluence.qps.nl/qinsy/en/world-geodetic-system-1984-wgs84-29855173.html

Information on vertical datums:
https://www.nrcan.gc.ca/earth-sciences/geomatics/geodetic-reference-systems/9054#_Toc372901506

One service known to work with gpsrinex output is CSRS-PPP at NRCAN:
https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php

Another service known to work with gpsrinex output is Trimble RTX
from Trimble.  They require dual frequency (L1 and L2) raw data:
https://trimblertx.com/Home.aspx

GAPS requires L2 P or L5 I+Q signals, and is not supported by gpsd:
http://gaps.gge.unb.ca/

OPUS requires L1/L2 frequency observation files, and has limited geographic
coverage:
https://www.ngs.noaa.gov/OPUS/

The curious can find the RINEX 3.04 format described here:
ftp://ftp.igs.org/pub/data/format/rinex304.pdf

// vim: set syntax=asciidoc: