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

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)

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)

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

pd.options.display.max_columns = 50
GNSSData.head(10)

nvalues = int(GNSSData.count()[0])
GNSSData.count()

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')

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)}

# Display GPS Heatmap from Disk
gpsheatmap = ImageDisp(filename=gpsmapfile)
gpsheatmap

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)

print('Done')