Get Lat and Long position
electronicDuck opened this issue · 7 comments
Hi there,
Thanks for the great software!
I’m very new to python so sorry if it’s a silly question.
I am hoping to use this package to monitor a base station. My thinking was use this package to send the RTK correction data to a F9P. Next was to read the corrected lat/long and compare where the receiver thinks it is to where it should be. Calculate the distance and then log the difference.
How do I read the position? Your rtk_example.py works perfectly, and I get the serial stream. If I turn on print full it’s all there. How do I save the lat and long as variables though?
I have tried a few things but have only managed to mangle the script so far and not achieve anything :P
Thank you
Hi @electronicDuck ,
So it's a relatively simple matter to output the lat/lon from the receiver. Establishing the deviation between these and a nominated set of fixed reference coordinates is potentially non-trivial but a simple method is to use the so-called "haversine" formula to calculate the spherical distance (i.e. distance on the earth's surface) between the receiver's coordinates and your fixed reference coordinates.
I've attached below a slightly modified version of the rtk_example.py script which includes a highlighted block of code to do this. As I say, this is not necessarily the most geospatially rigorous technique but it hopefully gives you some pointers.
As an alternative, you could monitor the hAcc (horizontal accuracy) attribute in messages like PUBX00 (NMEA) or NAV-PVT (UBX), which would give you an indirect indication of RTK positional accuracy.
Bear in mind you need an accurate fixed reference and a good external antenna with a clear view of the sky to get the best results.
Hope this helps - any queries, get back to me. I'll keep the issue open for a week or so.
"""
pygnssutils - rtk_example_deviation.py
*** FOR ILLUSTRATION ONLY - NOT FOR PRODUCTION USE ***
PLEASE RESPECT THE TERMS OF USE OF ANY NTRIP SERVER YOU
USE WITH THIS EXAMPLE - INAPPROPRIATE USE CAN RESULT IN
YOUR NTRIP USER ACCOUNT OR IP BEING TEMPORARILY BLOCKED.
This example illustrates how to use the UBXReader and
GNSSNTRIPClient classes to get RTCM3 RTK data from a
designated NTRIP server and apply it to an RTK-compatible
GNSS receiver (e.g. ZED-F9P) connected to a local serial port.
GNSSNTRIPClient receives RTK data from the NTRIP server
and outputs it to a message queue. A send thread reads data
from this message queue and sends it to the receiver. A read
thread reads and parses data from the receiver and prints it
to the terminal.
For brevity, the example will print out just the identities of
all incoming GNSS and NTRIP messages, but the full message can
be printed by setting the global PRINT_FULL variable to True.
If the receiver is a u-blox UBX receiver, it can be configured
to output RXM-RTCM messages which acknowledge receipt of
incoming RTK data and confirm whether or not it was used (i.e.
RTK correction applied).
NB: Some NTRIP servers may stop sending RTK data after a while
if they're not receiving legitimate NMEA GGA position updates
from the client.
Created on 5 Jun 2022
@author: semuadmin
"""
# pylint: disable=broad-except
from sys import platform
from io import BufferedReader
from threading import Thread, Lock, Event
from queue import Queue
from time import sleep
from serial import Serial
from pyubx2 import (
UBXReader,
NMEA_PROTOCOL,
UBX_PROTOCOL,
RTCM3_PROTOCOL,
protocol,
)
from pyrtcm import RTCM_MSGIDS
from pygnssutils import GNSSNTRIPClient, VERBOSITY_LOW, haversine
# NTRIP server parameters - AMEND AS REQUIRED:
# Ideally, mountpoint should be <30 km from location.
NTRIP_SERVER = "rtk2go.com"
NTRIP_PORT = 2101
MOUNTPOINT = "cbrookf9p"
NTRIP_USER = "anon"
NTRIP_PASSWORD = "password"
# Fixed reference coordinates - AMEND AS REQUIRED:
REFLAT = -33.4
REFLON = 138.2
REFALT = 124
REFSEP = 0
# Set to True to print entire GNSS/NTRIP message rather than just identity
PRINT_FULL = False
def read_gnss(stream, lock, stopevent):
"""
THREADED
Reads and parses incoming GNSS data from receiver.
"""
ubr = UBXReader(
BufferedReader(serial),
protfilter=(NMEA_PROTOCOL | UBX_PROTOCOL | RTCM3_PROTOCOL),
)
while not stopevent.is_set():
try:
if stream.in_waiting:
lock.acquire()
(raw_data, parsed_data) = ubr.read() # pylint: disable=unused-variable
lock.release()
if parsed_data:
# START OF DEVIATION LOGGING ILLUSTRATION
# if message contains lat/lon, compare it to reference coordinates
# uses haversine formula to estimate spherical distance ('deviation')
# between the two sets of coordinates
if hasattr(parsed_data, "lat") and hasattr(parsed_data, "lon"):
lat = parsed_data.lat
lon = parsed_data.lon
dev = haversine(lat, lon, REFLAT, REFLON) * 1000 # meters
print(
f"Receiver coordinates: {lat}, {lon}. Deviation from fixed ref: {dev:06f} m"
)
# END OF DEVIATION LOGGING ILLUSTRATION
idy = parsed_data.identity
# if it's an RXM-RTCM message, show which RTCM3 message
# it's acknowledging and whether it's been used or not.""
if idy == "RXM-RTCM":
nty = (
f" - {parsed_data.msgType} "
f"{'Used' if parsed_data.msgUsed > 0 else 'Not used'}"
)
else:
nty = ""
if PRINT_FULL:
print(parsed_data)
else:
print(f"GNSS>> {idy}{nty}")
except Exception as err:
print(f"Something went wrong in read thread {err}")
break
def send_gnss(stream, lock, stopevent, inqueue):
"""
THREADED
Reads RTCM3 data from message queue and sends it to receiver.
"""
while not stopevent.is_set():
try:
raw_data, parsed_data = inqueue.get()
if protocol(raw_data) == RTCM3_PROTOCOL:
if PRINT_FULL:
print(parsed_data)
else:
print(
f"NTRIP>> {parsed_data.identity} {RTCM_MSGIDS[parsed_data.identity]}"
)
lock.acquire()
stream.write(raw_data)
lock.release()
except Exception as err:
print(f"Something went wrong in send thread {err}")
break
def start_read_thread(stream, lock, stopevent):
"""
Start read thread.
"""
rth = Thread(
target=read_gnss,
args=(
stream,
lock,
stopevent,
),
daemon=True,
)
rth.start()
return rth
def start_send_thread(stream, lock, stopevent, msgqueue):
"""
Start send thread.
"""
kth = Thread(
target=send_gnss,
args=(
stream,
lock,
stopevent,
msgqueue,
),
daemon=True,
)
kth.start()
return kth
if __name__ == "__main__":
# pylint: disable=invalid-name
# GNSS receiver serial port parameters - AMEND AS REQUIRED:
if platform == "win32": # Windows
SERIAL_PORT = "COM13"
elif platform == "darwin": # MacOS
SERIAL_PORT = "/dev/tty.usbmodem1101"
else: # Linux
SERIAL_PORT = "/dev/ttyACM1"
BAUDRATE = 9600
TIMEOUT = 0.1
GGAMODE = 1 # use fixed reference position (0 = use live position)
GGAINT = 10 # -1 = do not send NMEA GGA sentences
serial_lock = Lock()
ntrip_queue = Queue()
stop = Event()
try:
print(f"Opening serial port {SERIAL_PORT} @ {BAUDRATE}...\n")
with Serial(SERIAL_PORT, BAUDRATE, timeout=TIMEOUT) as serial:
stop.clear()
print("Starting read thread...\n")
gnss_thread = start_read_thread(serial, serial_lock, stop)
print("Starting send thread...\n")
rtk_thread = start_send_thread(serial, serial_lock, stop, ntrip_queue)
print(f"Starting NTRIP client on {NTRIP_SERVER}:{NTRIP_PORT}...\n")
with GNSSNTRIPClient(None, verbosity=VERBOSITY_LOW) as gnc:
streaming = gnc.run(
server=NTRIP_SERVER,
port=NTRIP_PORT,
mountpoint=MOUNTPOINT,
user=NTRIP_USER,
password=NTRIP_PASSWORD,
reflat=REFLAT,
reflon=REFLON,
refalt=REFALT,
refsep=REFSEP,
ggamode=GGAMODE,
ggainterval=GGAINT,
output=ntrip_queue,
)
while streaming and not stop.is_set(): # run until user presses CTRL-C
sleep(1)
sleep(1)
except KeyboardInterrupt:
stop.set()
print("Terminated by user")
Thank you so much! Code works perfectly first time I ran it! Does exactly what I was hoping it would do :) Thank you for the change log, interesting to see which bits needed changing, has helped my understanding!
Only question I have is I am hoping the distance difference will be in the mm/cm range, is the Haversine method best for this or would a different technique be better?
Hi again @electronicDuck
Glad it helped, but to answer your question - no, the haversine formula used here is based on a rough approximation of the Earth as a perfect sphere, which is not exactly true at a local geographical level (though not far off - about 0.5% accurate). You could in theory apply local 'correction tables' to the formula to obtain slightly better accuracy, but my example was only offered by way of a simple illustration.
You may perhaps need to manage your expectations here. Whilst the F9P is capable of delivering cm level accuracy in ideal conditions, obtaining this level of accuracy is somewhat challenging in practice. Forget mm level accuracy - you're almost certainly not going to get that with non-professional equipment.
- Firstly, you'd need to know your basestation location to an equivalent level of accuracy, which generally requires specialist surveying techniques (no conventional map will give you this level of accuracy).
- It requires a good external antenna and ground plane with an unrestricted view of the sky. The chip antennas supplied with many hobbyist GNSS receivers or devkits are unlikely to deliver adequate results.
- It requires a reliable local DGPS (NTRIP) data source, ideally no more than 30km from the basestation and preferably closer (i.e. subject to the same local atmospheric conditions).
- Calculating the true surface distance between your receiver coordinates and the basestation coordinates - in so far as this is meaningful - would require an extremely accurate model of your local topography.
The following might be of interest if you're new to GNSS:
https://www.semuconsulting.com/gnsswiki/
I'll close this issue for now. I recommend you Google 'RTK' or 'DGPS', or consult the u-blox support forum, to obtain more detailed information on how best to utilise RTK techniques to improve GNSS accuracy.
By 'Cartesian position" I'm assuming you mean the Earth Centered Earth Fixed (ECEF) coordinates? These are indeed available from UBX messages e.g. NAV-POSECEF or NAV-HPPOSECEF (high precision), assuming the F9P has been configured to output these UBX message types - by default it only outputs a selection of NMEA messages. You could use u-Center or PyGPSClient to enable these messages.
But if you're measuring accuracy in terms of straight line deviation from a nominated 3D fixed reference point, you'd obviously need the accurate ECEF coordinates of your fixed reference point as well. Then it's just a case of applying a simple geometric calculation, but I'm not sure this would necessarily give you better results. It all depends on what 'accuracy' means in the context of your intended application.
From the description of your setup, I see no reason why it should not be possible to achieve "cm level" accuracy (i.e. 1-3 cm at best) from the F9P provided your DGPS source (NTRIP caster or other ground-based reference) is in close enough proximity. The most significant error in modern GNSS measurements generally comes from ionospheric effects, and this is essentially what DGPS techniques are compensating for.
"mm level" accuracy relative to a terrestrial geodetic datum (e.g. WGS84) would involve taking into account such factors as tectonic plate shift, local gravitational anomalies, etc. etc. It's simply not realistic using mainstream GNSS technology.
As I mentioned above, if you want an alternative indication of nominal horizontal accuracy, you could refer to the hAcc attribute in messages like PUBX00 (NMEA) or NAV-PVT (UBX) - again, you'd need to configure your F9P to output these message types, and add the following lines to the modified rtk_example.py script:
if hasattr(parsed_data, "hAcc"):
print(f"Horizontal accuracy: {parsed_data.hAcc}")
Hope this helps.