import pandas as pd
import numpy as np
from datetime import datetime, timedelta
from IPython.display import Image as ImageDisp
from pandas import DataFrame
import string
import os
import matplotlib.pyplot as plt
%pylab inline --no-import-all
Populating the interactive namespace from numpy and matplotlib
def nmeatime(date):
return datetime.strptime(date, "%H%M%S.%f")
fto = './DataLogs/Messfahrt_2_start64.log' # File to Open
data = pd.read_table(fto, sep=',', header=None, index_col=1, parse_dates=True, date_parser=nmeatime, comment='*')
GSTData = data[1::2].drop(0,1)
1 2 3 4 5 6 7 8
| | | | | | | |
$--GST,hhmmss.ss,x.x,x.x,x.x,x.x,x.x,x.x,x.x*hh<CR><LF>
This message is used to support Receiver Autonomous Integrity Monitoring (RAIM). Pseudorange measurement error statistics can be translated in the position domain in order to give statistical measures of the quality of the position solution. If only GPS, GLONASS, etc. is used for the reported position solution, the talker ID is GP, GL, etc., and the error data pertains to the individual system. If satellites from multiple systems are used to obtain the reported position solution, the talker ID is GN and the errors pertain to the combined solution.
GSTData.rename(columns={1: 'UTC Time', 2: 'RMS', 3: 'SigmaMajorAxis', 4: 'SigmaMinorAxis', 5: 'OrientMajorAxis', 6: 'SigmaLat', 7: 'SigmaLon', 8: 'SigmaAlt'}, inplace=True)
GSTData=GSTData.drop([9, 10, 11, 12, 13, 14],1)
GSTData.head(10)
RMS | SigmaMajorAxis | SigmaMinorAxis | OrientMajorAxis | SigmaLat | SigmaLon | SigmaAlt | |
---|---|---|---|---|---|---|---|
1 | |||||||
1900-01-01 08:09:06.400000 | 0.786 | 0.141 | 0.022 | 159.752 | 0.133 | 0.053 | 0.772 |
1900-01-01 08:09:06.500000 | 0.788 | 0.142 | 0.022 | 159.752 | 0.133 | 0.053 | 0.775 |
1900-01-01 08:09:06.600000 | 0.766 | 0.138 | 0.021 | 159.752 | 0.130 | 0.052 | 0.753 |
1900-01-01 08:09:06.700000 | 0.765 | 0.138 | 0.021 | 159.810 | 0.129 | 0.051 | 0.753 |
1900-01-01 08:09:06.800000 | 0.763 | 0.137 | 0.021 | 159.761 | 0.129 | 0.051 | 0.751 |
1900-01-01 08:09:06.900000 | 0.766 | 0.138 | 0.021 | 159.761 | 0.129 | 0.051 | 0.753 |
1900-01-01 08:09:07 | 0.768 | 0.138 | 0.021 | 159.761 | 0.130 | 0.052 | 0.755 |
1900-01-01 08:09:07.100000 | 0.771 | 0.139 | 0.021 | 159.761 | 0.130 | 0.052 | 0.758 |
1900-01-01 08:09:07.200000 | 0.773 | 0.139 | 0.021 | 159.761 | 0.131 | 0.052 | 0.760 |
1900-01-01 08:09:07.300000 | 0.775 | 0.139 | 0.021 | 159.761 | 0.131 | 0.052 | 0.762 |
11
1 2 3 4 5 6 7 8 9 10 | 12 13 14 15
| | | | | | | | | | | | | | |
$--GGA,hhmmss.ss,llll.ll,a,yyyyy.yy,a,x,xx,x.x,x.x,M,x.x,M,x.x,xxxx*hh<CR><LF>
Field Number:
0 - fix not available
1 - GPS fix
2 - Differential GPS fix
GGAData = data[::2].drop(0,1)
GGAData.rename(columns={1: 'UTC Time', 2: 'Lat', 3: 'NorS', 4: 'Lon', 5: 'EorW', 6: 'Fix', 7: 'NumSats', 8: 'HDOP', 9: 'Alt', 10: 'Unit', 11: 'GeoidSeparation', 12: 'UnitGeoSep', 13: 'DGPSage', 14: 'DGPSStation'}, inplace=True)
GGAData['LatDD'] = (GGAData.Lat/100).fillna(0).astype(int)
GGAData['LatDD'] = GGAData.LatDD + (GGAData.Lat - 100.0*GGAData.LatDD)/60.0
GGAData['LonDD'] = (GGAData.Lon/100).fillna(0).astype(int)
GGAData['LonDD'] = GGAData.LonDD + (GGAData.Lon - 100.0*GGAData.LonDD)/60.0
GNSSData = pd.merge(GSTData, GGAData, left_index=True, right_index=True, how='outer')
GNSSData.index.names = ['UTC Time']
GNSSData.index = GNSSData.index.map(lambda t: t.replace(year=2014, month=4, day=23)) - pd.offsets.Hour(-2)
GNSSData.index
<class 'pandas.tseries.index.DatetimeIndex'> [2014-04-23 10:09:06.400000, ..., 2014-04-23 10:21:20.400000] Length: 7341, Freq: None, Timezone: None
pd.options.display.max_columns = 50
GNSSData.head(10)
RMS | SigmaMajorAxis | SigmaMinorAxis | OrientMajorAxis | SigmaLat | SigmaLon | SigmaAlt | Lat | NorS | Lon | EorW | Fix | NumSats | HDOP | Alt | Unit | GeoidSeparation | UnitGeoSep | DGPSage | DGPSStation | LatDD | LonDD | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
2014-04-23 10:09:06.400000 | 0.786 | 0.141 | 0.022 | 159.752 | 0.133 | 0.053 | 0.772 | 5100.753576 | N | 1347.488194 | E | 5 | 15 | 0.8 | 120.3053 | M | 43.9749 | M | 1.4 | 417 | 51.01256 | 13.79147 |
2014-04-23 10:09:06.500000 | 0.788 | 0.142 | 0.022 | 159.752 | 0.133 | 0.053 | 0.775 | 5100.753578 | N | 1347.488196 | E | 5 | 15 | 0.8 | 120.3183 | M | 43.9749 | M | 1.5 | 417 | 51.01256 | 13.79147 |
2014-04-23 10:09:06.600000 | 0.766 | 0.138 | 0.021 | 159.752 | 0.130 | 0.052 | 0.753 | 5100.753578 | N | 1347.488195 | E | 5 | 15 | 0.8 | 120.3022 | M | 43.9749 | M | 0.6 | 417 | 51.01256 | 13.79147 |
2014-04-23 10:09:06.700000 | 0.765 | 0.138 | 0.021 | 159.810 | 0.129 | 0.051 | 0.753 | 5100.753615 | N | 1347.488176 | E | 5 | 15 | 0.8 | 119.9284 | M | 43.9749 | M | 0.7 | 417 | 51.01256 | 13.79147 |
2014-04-23 10:09:06.800000 | 0.763 | 0.137 | 0.021 | 159.761 | 0.129 | 0.051 | 0.751 | 5100.753594 | N | 1347.488189 | E | 5 | 15 | 0.8 | 120.2211 | M | 43.9749 | M | 0.8 | 417 | 51.01256 | 13.79147 |
2014-04-23 10:09:06.900000 | 0.766 | 0.138 | 0.021 | 159.761 | 0.129 | 0.051 | 0.753 | 5100.753594 | N | 1347.488187 | E | 5 | 15 | 0.8 | 120.2124 | M | 43.9749 | M | 0.9 | 417 | 51.01256 | 13.79147 |
2014-04-23 10:09:07 | 0.768 | 0.138 | 0.021 | 159.761 | 0.130 | 0.052 | 0.755 | 5100.753594 | N | 1347.488186 | E | 5 | 15 | 0.8 | 120.2028 | M | 43.9749 | M | 1.0 | 417 | 51.01256 | 13.79147 |
2014-04-23 10:09:07.100000 | 0.771 | 0.139 | 0.021 | 159.761 | 0.130 | 0.052 | 0.758 | 5100.753594 | N | 1347.488187 | E | 5 | 15 | 0.8 | 120.2140 | M | 43.9749 | M | 1.1 | 417 | 51.01256 | 13.79147 |
2014-04-23 10:09:07.200000 | 0.773 | 0.139 | 0.021 | 159.761 | 0.131 | 0.052 | 0.760 | 5100.753595 | N | 1347.488187 | E | 5 | 15 | 0.8 | 120.2164 | M | 43.9749 | M | 1.2 | 417 | 51.01256 | 13.79147 |
2014-04-23 10:09:07.300000 | 0.775 | 0.139 | 0.021 | 159.761 | 0.131 | 0.052 | 0.762 | 5100.753594 | N | 1347.488187 | E | 5 | 15 | 0.8 | 120.2140 | M | 43.9749 | M | 1.3 | 417 | 51.01256 | 13.79147 |
nvalues = int(GNSSData.count()[0])
GNSSData.count()
RMS 7272 SigmaMajorAxis 7272 SigmaMinorAxis 7272 OrientMajorAxis 7272 SigmaLat 7272 SigmaLon 7272 SigmaAlt 7272 Lat 7272 NorS 7272 Lon 7272 EorW 7272 Fix 7341 NumSats 7341 HDOP 7272 Alt 7272 Unit 7272 GeoidSeparation 7272 UnitGeoSep 7272 DGPSage 6301 DGPSStation 6301 LatDD 7272 LonDD 7272 dtype: int64
fw = int(nvalues/6000.0)
if fw < 16:
fw=16
elif fw > 3000:
fw=3000
GNSSData['tvalue'] = GNSSData.index
GNSSData['delta'] = (GNSSData['tvalue'] - GNSSData['tvalue'].shift()).fillna(0)
dt_s = GNSSData['delta'].values.astype('int')/1e9
GNSSData['fa']=np.divide(1.0,dt_s)
print('Abtastrate der Messung (im Mittel): %.1fHz' % np.divide(1.0,dt_s.mean()))
plt.figure(figsize=(fw, 4))
GNSSData['fa'].plot()
plt.ylabel('$Hz$')
plt.title('Abtastrate')
#plt.axhline(50, color='k', alpha=0.2)
#plt.savefig(outdir+fto[16:-8]+'fa.png', dpi=150, bbox_inches='tight')
Abtastrate der Messung (im Mittel): 10.0Hz
<matplotlib.text.Text at 0x104731310>
gpslogfile = fto[:-4]+'-LatLon.log'
gpsmapfile = fto[:-4]+'-GPSHeatmap.png'
gpsmapwidth= 800 # px
gpsmapmargin=100 # px
gpsmapbase = 'http://b.tile.stamen.com/toner-lite'
# write lat/lon text file to read in to heatmap
GNSSData[['LatDD','LonDD']][::10].to_csv(gpslogfile, sep=' ', index=False, header=False)
# create Heatmap of Lat/Lon
# -M =Value Color in HSV Hex with Direction DHHSSVVAA (=Direction, Hue, Suration, Value, Alpha)
%run heatmap.py {'-p %s -o %s --width %i -M 018FFFFFF --margin %i --osm --osm_base %s' % (gpslogfile, gpsmapfile, gpsmapwidth, gpsmapmargin, gpsmapbase)}
WARNING: Using /var/folders/yz/7x09z5gd5dg_4q1mkf9clsfr0000gn/T/ to cache maptiles. Retrieving 24 tiles... ... done.
# Display GPS Heatmap from Disk
gpsheatmap = ImageDisp(filename=gpsmapfile)
gpsheatmap
Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under CC BY SA.
plt.figure(figsize=(fw, 4))
la=GNSSData['RMS'].plot(alpha=0.6)
la=GNSSData['DGPSage'].plot(alpha=0.6)
plt.legend(loc=2)
ra=GNSSData.NumSats.plot(secondary_y=True, alpha=0.4)
la.set_ylabel('RMS ($m$) und DGPS Age')
#la.set_ylim([-6,6])
ra.right_ax.set_ylabel('Satellites Used')
#ra.set_ylim([0, 8])
plt.legend(loc=1)
plt.savefig(fto[:-4]+'-RMS.png', dpi=150, bbox_inches='tight')
whratio = np.cos(np.mean(GNSSData.LatDD*np.pi/180.0))
fh = 14.0
fig = plt.figure(figsize=(fh, fh*whratio))
# Measurements
every = 50 #
plt.scatter(GNSSData.LonDD[::every], GNSSData.LatDD[::every], s=100, label='GNSS measurement (every %sst)' % every,\
c=GNSSData.RMS[::every], cmap='autumn_r')
cbar=plt.colorbar()
cbar.ax.set_ylabel(u'RMS', rotation=270)
cbar.ax.set_xlabel(u'$m$')
plt.xlabel('Longitude [$^\circ$]')
plt.ylabel('Latitude [$^\circ$]')
plt.legend(loc='best')
plt.axis('equal')
#plt.tight_layout()
# xticks
locs,labels = plt.xticks()
plt.xticks(locs, map(lambda x: "%.3f" % x, locs));
# ytikcs
locs,labels = plt.yticks()
plt.xticks(rotation=35)
plt.yticks(locs, map(lambda x: "%.3f" % x, locs));
plt.savefig(fto[:-4]+'-RMS-LatLon.png', dpi=150, bbox_inches='tight')
plt.figure(figsize=(fw, 4))
la=GNSSData['Alt'].plot(alpha=0.6)
plt.xlabel('Time')
plt.legend(loc=2)
ra=GNSSData.SigmaAlt.plot(secondary_y=True, alpha=0.4)
la.set_ylabel('Altitude ($m$)')
#la.set_ylim([-6,6])
ra.right_ax.set_ylabel('$\sigma_{Altitude}$')
#ra.set_ylim([0, 8])
plt.legend(loc=1)
<matplotlib.legend.Legend at 0x106311790>
print('Done')
Done