)
Trimble NetR9 Receiver - GSOF Messages
General Serial Output Format (GSOF) messages are a Trimble proprietary format and can be used to send information such as position and status to a third-party device.
More detailed information regarding the structure of this binary message format can be found at Trimble’s GSOF Specification Page.
To enable a GSOF stream on a Trimble NetR9 receiver:
- Go to ’I/O Configuration’ on the web interface.
- Either select an unused serial port or ’add TCP/IP or UDP Port.
- Select the ’GSOF’ message type from the pull down menu.
- If using ’add TCP/IP or UDP Port’, specify the port number that you’d like the stream to broadcast from.
- If using a serial port, configure the baud, parity and flow to match your third-pary device.
- Select which messages you would like the receiver to output.
- Click ’OK’ to enable the stream.
See below for an example python scripts ( getGsofStream.py and utilities.py) connect to a TCP socket and parses GSOF message.
Sample output:
Reference not specified, using first epoch: X =-1282458.7777 Y=-4718377.8575 Z=4084176.9870
WN - GPS WEEK
SOW - GPS time
FLAG1 FLAG2
SATS - # SVS
X[m] Y[m] Z[m]
eX[m] eY[m] eZ[m] - Position VCV Position RMS
E[m] N[m] U[m]
# WN SOW FLAG1 FLAG2 SATS X[m] Y[m] Z[m] eX[m] eY[m] eZ[m] E[m] N[m] U[m]
2035 319837.000 191 35 16 -1282458.7759 -4718377.8519 4084176.9822 0.0144 0.0232 0.0174 0.0002 0.0001 -0.0075
2035 319838.000 191 35 16 -1282458.7719 -4718377.8502 4084176.9807 0.0146 0.0234 0.0176 0.0037 0.0007 -0.0106
2035 319839.000 191 35 16 -1282458.7801 -4718377.8574 4084176.9830 0.0148 0.0237 0.0178 -0.0024 -0.0034 -0.0021
############### getGsofStream.py ###################
from struct import unpack
from utilities import *
import socket
import math
import argparse
__author__ = "Henry T. Berglund"
class Gsof(object):
""" A class to parse GSOF messages from a TCP stream. """
def __init__(self, sock=None):
if sock is None:
self.sock = socket.socket(
socket.AF_INET, socket.SOCK_STREAM)
else:
self.sock = sock
self.msg_dict = {}
self.msg_bytes = None
self.checksum = None
self.rec_dict = {}
def connect(self, host, port):
self.sock.connect((host, port))
def get_message_header(self):
data = self.sock.recv(7)
msg_field_names = (’STX’, ’STATUS’, ’TYPE’, ’LENGTH’,
’T_NUM’, ’PAGE_INDEX’, ’MAX_PAGE_INDEX’)
self.msg_dict = dict(zip(msg_field_names, unpack(’>7B’, data)))
self.msg_bytes = self.sock.recv(self.msg_dict[’LENGTH’] - 3)
(checksum, etx) = unpack(’>2B’, self.sock.recv(2))
def checksum256(st):
"""Calculate checksum"""
return reduce(lambda x, y: x+y, map(ord, st)) % 256
if checksum-checksum256(self.msg_bytes+data[1:]) == 0:
self.checksum = True
else:
self.checksum = False
def get_records(self):
while len(self.msg_bytes) > 0:
# READ THE FIRST TWO BYTES FROM RECORD HEADER
record_type, record_length = unpack(’>2B’, self.msg_bytes[0:2])
self.msg_bytes = self.msg_bytes[2:]
self.select_record(record_type, record_length)
def select_record(self, record_type, record_length):
if record_type == 1:
rec_field_names = (’GPS_TIME’, ’GPS_WEEK’, ’SVN_NUM’,
’FLAG_1’, ’FLAG_2’, ’INIT_NUM’)
rec_values = unpack(’>LH4B’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 2:
rec_field_names = (’LATITUDE’, ’LONGITUDE’, ’HEIGHT’)
rec_values = unpack(’>3d’, self.msg_bytes[0:record_length])
rec_values = list(rec_values)
rec_values = (math.degrees(rec_values[0]), math.degrees(rec_values[1]), rec_values[2])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 3:
rec_field_names = (’X_POS’, ’Y_POS’, ’Z_POS’)
rec_values = unpack(’>3d’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 4:
rec_field_names = (’LOCAL_DATUM_ID’, ’LOCAL_DATUM_LAT’,
’LOCAL_DATUM_LON’, ’LOCAL_DATUM_HEIGHT’, ’OPRT’)
rec_values = unpack(’>8s3dB’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 5:
rec_field_names = (’LOCAL_DATUM_ID’, ’LOCAL_ZONE_ID’,
’LOCAL_ZONE_NORTH’, ’LOCAL_ZONE_EAST’, ’LOCAL_DATUM_HEIGHT’)
rec_values = unpack(’>2s3d’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 6:
rec_field_names = (’DELTA_X’, ’DELTA_Y’, ’DELTA_Z’)
rec_values = unpack(’>3d’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 7:
rec_field_names = (’DELTA_EAST’, ’DELTA_NORTH’, ’DELTA_UP’)
rec_values = unpack(’>3d’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 8:
rec_field_names = (’VEL_FLAG’, ’VELOCITY’, ’HEADING’, ’VERT_VELOCITY’)
rec_values = unpack(’>B3f’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 9:
rec_field_names = (’PDOP’, ’HDOP’, ’VDOP’, ’TDOP’)
rec_values = unpack(’>4f’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 10:
rec_field_names = (’CLOCK_FLAG’, ’CLOCK_OFFSET’, ’FREQ_OFFSET’)
rec_values = unpack(’>B2d’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 11:
rec_field_names = (’POSITION_RMS_VCV’, ’VCV_XX’, ’VCV_XY’, ’VCV_XZ’,
’VCV_YY’, ’VCV_YZ’, ’VCV_ZZ’, ’UNIT_VAR_VCV’, ’NUM_EPOCHS_VCV’)
rec_values = unpack(’>8fh’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 12:
rec_field_names = (’POSITION_RMS_SIG’, ’SIG_EAST’, ’SIG_NORT’, ’COVAR_EN’, ’SIG_UP’,
’SEMI_MAJOR’, ’SEMI_MINOR’, ’ORIENTATION’, ’UNIT_VAR_SIG’,
’NUM_EPOCHS_SIG’)
rec_values = unpack(’>9fh’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 15:
rec_field_names = ’SERIAL_NUM’
rec_values = unpack(’>l’, self.msg_bytes[0:record_length])
self.rec_dict.update({rec_field_names: rec_values[0]})
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 16:
rec_field_names = (’GPS_MS_OF_WEEK’, ’CT_GPS_WEEK’, ’UTC_OFFSET’, ’CT_FLAGS’)
rec_values = unpack(’>l2hB’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 26:
rec_field_names = (’UTC_MS_OF_WEEK’, ’UTC_GPS_WEEK’, ’UTC_SVS_NUM’, ’UTC_FLAG_1’, ’UTC_FLAG_2’,
’UTC_INIT_NUM’)
rec_values = unpack(’>lh4B’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 34:
NUM_OF_SVS = unpack(’>B’, self.msg_bytes[0])
self.msg_bytes = self.msg_bytes[1:]
rec_field_names = (’PRN’, ’SV_SYSTEM’, ’SV_FLAG1’, ’SV_FLAG2’, ’ELEVATION’, ’AZIMUTH’,
’SNR_L1’, ’SNR_L2’, ’SNR_L5’)
for field in xrange(len(rec_field_names)):
self.rec_dict[rec_field_names[field]] = []
for sat in xrange(NUM_OF_SVS[0]):
rec_values = unpack(’>5Bh3B’, self.msg_bytes[0:10])
self.msg_bytes = self.msg_bytes[10:]
for num in xrange(len(rec_field_names)):
self.rec_dict[rec_field_names[num]].append(rec_values[num])
elif record_type == 37:
rec_field_names = (’BATT_CAPACITY’, ’REMAINING_MEM’)
rec_values = unpack(’>hd’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
elif record_type == 41:
rec_field_names = (’GPS_TIME_REF’, ’GPS_WEEK_REF’, ’LAT_REF’, ’LONG_REF’,
’HEIGHT_REF’, "QI")
rec_values = unpack(’>LH3dB’, self.msg_bytes[0:record_length])
self.rec_dict.update(dict(zip(rec_field_names, rec_values)))
self.msg_bytes = self.msg_bytes[record_length:]
else:
"""Unknown record type? Skip it for now!"""
#print record_type
self.msg_bytes = self.msg_bytes[record_length:]
def main():
parser = argparse.ArgumentParser(prog=’getGsofStream’, usage=’%(prog)s -i IP_ADDRESS -p PORT -r X Y Z ’)
parser.add_argument(’-i’, ’--IP’, action="store", default="69.44.87.212", type=str)
parser.add_argument(’-p’, ’--PORT’, action="store", default=9999, type=int)
parser.add_argument(’-r’, ’--ref’, nargs=3, action=’append’)
results = parser.parse_args()
# OPEN GSOF STREAM
s = Gsof()
s.connect(results.IP, results.PORT)
try:
# Use user provided reference coordinates if given
ref_x, ref_y, ref_z = results.ref[0]
print "# Reference specified, using X =%.3f Y=%.3f Z=%.3f" % (ref_x, ref_y, ref_z)
except:
# Get reference coordinates from first epoch if not provided by user
s.get_message_header()
s.get_records()
ref_x, ref_y, ref_z = s.rec_dict[’X_POS’], s.rec_dict[’Y_POS’], s.rec_dict[’Z_POS’]
print "# Reference not specified, using first epoch: X =%.4f Y=%.4f Z=%.4f" % (ref_x, ref_y, ref_z)
print ""
print "WN - GPS WEEK"
print "SOW - GPS time"
print "FLAG1 FLAG2"
print "SATS - # SVS"
print "X[m] Y[m] Z[m]"
print "eX[m] eY[m] eZ[m] - Position VCV Position RMS"
print "E[m] N[m] U[m]"
print ""
print "# WN SOW FLAG1 FLAG2 SATS X[m] Y[m] Z[m] eX[m] eY[m] eZ[m] E[m] N[m] U[m]"
while 1:
# READ GSOF STREAM
s.get_message_header()
s.get_records()
# PRINT GSOF STREAM
x = s.rec_dict[’X_POS’]
y = s.rec_dict[’Y_POS’]
z = s.rec_dict[’Z_POS’]
neu = xyz2neu(ref_x, ref_y, ref_z, x, y, z)
output = "%04d %.3f %3d %3d %2d %14.4f %14.4f %14.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f" % (
s.rec_dict[’GPS_WEEK’],
s.rec_dict[’GPS_TIME’]/1000.0,
s.rec_dict[’FLAG_1’],
s.rec_dict[’FLAG_2’],
s.rec_dict[’SVN_NUM’],
s.rec_dict[’X_POS’],
s.rec_dict[’Y_POS’],
s.rec_dict[’Z_POS’],
math.sqrt(s.rec_dict[’VCV_XX’]),
math.sqrt(s.rec_dict[’VCV_YY’]),
math.sqrt(s.rec_dict[’VCV_ZZ’]),
neu[1],
neu[0],
neu[2],
)
print(output)
if __name__ == ’__main__’:
try:
main()
except KeyboardInterrupt:
print "\n#: Stopping, you hit ctl+C. "
#traceback.print_exc()
###### utilities.py #################
__author__ = ’berglund’
"""
Utilities for various geodetic conversions
"""
from math import sin, cos, sqrt, atan2, degrees, radians
def llh2xyz(lat, lon, h, a=6378137.0, finv=298.257223563, input_degrees=True):
"""
Converts ellipsoidal coordinates to cartesian.
:param lat:
:param lon:
:param h:
:return: A tuple containing cartesian coordinates in ECEF [X, Y, Z] meters.
"""
# convert degrees to radians if necessary
if input_degrees:
lat = radians(lat)
lon = radians(lon)
# Default (a=6378137.0, finv=298.257223563) is for WGS84
f = 1/finv
b = a*(1 - f)
e2 = 1 - (1 - f)**2
# Compute the Cartesian coordinates
v = a/sqrt(1-e2*sin(lat)*sin(lat))
x = (v+h)*cos(lat)*cos(lon)
y = (v+h)*cos(lat)*sin(lon)
z = (v*(1-e2)+h)*sin(lat)
return x, y, z
def xyz2llh(x, y, z, a=6378137.0, finv=298.257223563, output_degrees=True):
"""
Converts cartesian coordinates to ellipsoidal. Uses an iterative algorithm.
:param x: x coordinate in meters
:param y: y coordinate in meters
:param z: z coordinate in meters
:param a: major semi-axis (Default: WGS84)
:param finv: inverse of flattening (Default: WGS84)
:param output_degrees: A boolean to select output in degrees (default is True)
:return: A tuple containing geodetic coordinates in WGS84 (lat, lon, h).
WGS84: a=6378137.0, finv=298.257223563
"""
# Default (a=6378137.0, finv=298.257223563) is for WGS84
f = 1/finv
b = a*(1 - f)
e2 = 1 - (1 - f)**2
# Latitude and height convergence criteria
elat = 1.0e-12
eht = 1.e-5
# Initial values for iteration
p = sqrt(x*x + y*y)
lat = atan2(z, p*(1-e2))
h = 0
dh = 1
dlat = 1
# Iterate until lat & h converge to elat & eht
while dlat > elat or dh > eht:
lat0 = lat
h0 = h
v = a/sqrt(1-e2*sin(lat)*sin(lat))
h = p/cos(lat)-v
lat = atan2(z, p*(1-e2*v/(v+h)))
dlat = abs(lat-lat0)
dh = abs(h-h0)
lon = atan2(y, x)
if output_degrees:
# output degrees
geodetic_coordinates = (degrees(lat), degrees(lon), h)
else:
# output radians
geodetic_coordinates = (lat, lon, h)
return geodetic_coordinates
def xyz2neu(X, Y, Z, x, y, z):
"""
Rotates the <(x-X),(y-Y),(z-Z)> vector into a topocentric coordinate system
tangential to the reference point XYZ.
Assumes coordinates are given in meters.
:param X: Reference Coordinate X (meters)
:param Y: Reference Coordinate Y (meters)
:param Z: Reference Coordinate Z (meters)
:param x: Observed location x (meters)
:param y: Observed location y (meters)
:param z: Observed location z (meters)
:return: A tuple containing topocentric coordinates in n, e, u (meters)
WGS84: a=6378137.0, finv=298.257223563
"""
lat, lon, h = xyz2llh(X, Y, Z, a=6378137.0, finv=298.257223563, output_degrees=False)
e = -sin(lon) * (x - X) + cos(lon) * (y - Y)
n = -sin(lat) * cos(lon) * (x - X) - sin(lat) * sin(lon) * (y - Y) + cos(lat) * (z - Z)
u = cos(lat) * cos(lon) * (x - X) + cos(lat) * sin(lon) * (y - Y) + sin(lat) * (z - Z)
return n, e, u
def xyz_to_globk_neu(X, Y, Z, x, y, z):
"""
Finds the arc-length between two LLH points in dN and dE.
The parameter dU is defined as the difference between the observed and reference ellipsoidal height.
Assumes coordinates are given in meters.
:param X: Reference Coordinate X (meters)
:param Y: Reference Coordinate Y (meters)
:param Z: Reference Coordinate Z (meters)
:param x: Observed location x (meters)
:param y: Observed location y (meters)
:param z: Observed location z (meters)
:return: A tuple containing n, e, u (meters)
"""
def sub(a, b):
"""
Subtracts elements in two lists
"""
c = [i - j for i, j in zip(a, b)]
return c
def mul(a, b):
"""
Multiplies elements in two lists
"""
c = [i * j for i, j in zip(a, b)]
return c
P1 = xyz2llh(X, Y, Z, a=6378137.0, finv=298.257223563, output_degrees=False)
P2 = xyz2llh(x, y, z, a=6378137.0, finv=298.257223563, output_degrees=False)
# Subtract the reference coordinate from the observed coordinate
P_diff = sub(P2, P1)
# Convert the dLLH to dNEU
# Assumes a sphere
rot = [6378137.0, 6378137.0*cos(P1[0]), 1]
NEU = mul(P_diff, rot)
return tuple(NEU)
def covrot(X, Y, Z, cxx, cyy, czz, cxy, cxz, cyz):
"""
Rotates the sigmas and the correlations from xyz to enu
:param cxx:
:param cyy:
:param czz:
:param cxy:
:param cxz:
:param cyz:
:return:
"""
def matmult(a, b):
zip_b = zip(*b)
return [[sum(ele_a*ele_b for ele_a, ele_b in zip(row_a, col_b)) for col_b in zip_b] for row_a in a]
lat, lon, h = xyz2llh(X, Y, Z, a=6378137.0, finv=298.257223563, output_degrees=False)
rxx = -sin(lat)*cos(lon)
rxy = -sin(lat)*sin(lon)
rxz = cos(lat)
ryx = -sin(lon)
ryy = cos(lon)
ryz = 0
rzx = cos(lat)*cos(lon)
rzy = cos(lat)*sin(lon)
rzz = sin(lat)
rotation_matrix = [[rxx, rxy, rxz],
[ryx, ryy, ryz],
[rzx, rzy, rzz]]
rotation_transpose = [[rxx, ryx, rzx],
[rxy, ryy, rzy],
[rxz, ryz, rzz]]
covariance_matrix = [[cxx**2, cxy*cxx*cyy, cxz*cxx*czz],
[cxy*cxx*cyy, cyy**2, cyz*cyy*czz],
[cxz*cxx*czz, cyz*cyy*czz, czz**2]]
tmp = matmult(rotation_matrix, covariance_matrix)
result = matmult(tmp, rotation_transpose)
cxx = result[0][0]**0.5
cyy = result[1][1]**0.5
czz = result[2][2]**0.5
cxy = result[0][1]/(cxx*cyy)
cxz = result[0][2]/(cxx*czz)
cyz = result[1][2]/(cyy*czz)
return cxx, cyy, czz, cxy, cxz, cyz
Code disclaimer information
This document contains programming examples. UNAVCO grants you a nonexclusive copyright license to use all programming code examples from which you can generate similar function tailored to your own specific needs. All sample code is provided by UNAVCO for illustrative purposes only. These examples have not been thoroughly tested under all conditions. All programs contained herein are provided to you "AS IS" without any warranties of any kind.