#!/usr/bin/python # gps_hill_speed_plot.py # --copyright-- Copyright 2016 (C) Tranzoa, Co. All rights reserved. Warranty: You're free and on your own here. This code is not necessarily up-to-date or of public quality. # --url-- http://www.tranzoa.net/tzpython/ # --email-- pycode is the name to send to. tranzoa.com is the place to send to. # --bodstamps-- # December 20, 2016 bar # December 21, 2016 bar better labels # December 22, 2016 bar read/write speeds' json/pickle file # save image file # try --maximize (bad idea) # December 24, 2016 bar --exclude # multi-page pdf # better auto tick labelling # base Y at zero if all dots are above zero # better handle single day's track # December 25, 2016 bar don't combine multiple tracks per day. Per-Day -> Per-Track # November 7, 2017 bar maxint->maxsize # --eodstamps-- ## \file # \namespace tzpython.gps_hill_speed_plot # # # TODO: # GMT or local. Find the local time at the lat/lon? # Mouse-hover over a dot should show information about the dot - GPS track file name, for instance. # Be more robust to inputting real tracks rather than hikified, sparsified tracks. # # """ Plot point to point-ish speed as a function of when the points were and how steep the altitude change is between the points. """ import copy import datetime import fnmatch import json import math import os import sys import time import matplotlib.dates as mdates import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages import latlon import tz_gps import tz_parse_time import tzlib A_DAY = 24 * 3600 def point_count(speeds) : return(sum([ len(sa) for sa in speeds.values() ] + [ 0 ])) def on_day(t) : """ Return the day this time is in. """ return(int(t / A_DAY)) # GMT isn't right to use, really, maybe, though because we use only the first point in the track, maybe this is better than the localtime kludge tm = time.localtime(t) return(tm.tm_yday + (tm.tm_year * 366)) def plot_years(plt, sa, start_when, end_when, lbl) : ax = plt.gca() lct = mdates.AutoDateLocator() ax.xaxis.set_major_locator(lct) end_when += A_DAY tks = [] if (end_when - start_when) / A_DAY >= 365 * 1.5 : ll = 4 # ax.xaxis.set_major_formatter(mdates.AutoDateFormatter(lct)) ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y')) ax.xaxis.set_minor_locator(mdates.MonthLocator()) fy = time.localtime(start_when).tm_year ty = time.localtime(end_when ).tm_year + 1 while fy < ty : t = tz_parse_time.make_utime(1, 1, fy, 12, 0, 0) if start_when < t < end_when : # print "@@@@", fy, ty, t tks.append(datetime.datetime.fromtimestamp(t)) fy += 1 pass elif (end_when - start_when) / A_DAY >= 100 : tks = [] ll = 10 ax.xaxis.set_major_formatter(mdates.DateFormatter('%m/%d/%Y')) ax.xaxis.set_minor_locator(mdates.DayLocator()) etm = time.localtime(end_when) tm = time.localtime(start_when) y = tm.tm_year m = tm.tm_mon while True : if (y > etm.tm_year) or ((y >= etm.tm_year) and (m > etm.tm_mon)) : break t = tz_parse_time.make_utime(m, 1, y, 12, 0, 0) if start_when < t < end_when : tks.append(datetime.datetime.fromtimestamp(t)) m += 1 if m > 12 : y += 1 m = 1 pass pass else : tks = [] ll = 10 ax.xaxis.set_major_formatter(mdates.DateFormatter('%m/%d/%Y')) ax.xaxis.set_minor_locator(mdates.DayLocator()) tm = time.localtime(start_when) t = start_when - (tm.tm_wday * A_DAY) + (12 * 3600) while t < end_when : if start_when < t < end_when : tks.append(datetime.datetime.fromtimestamp(t)) t += 7 * A_DAY pass ax.xaxis.set_ticks(tks) ad = 'left' if len(tks) * ll > 48 : plt.gcf().autofmt_xdate() # note: this must be done here. Not earlier. Not later. ad = 'right' for tik in ax.xaxis.get_major_ticks() : tik.label1.set_horizontalalignment(ad) # # for when the cursor is hovering over the graph # ax.format_xdata = mdates.DateFormatter('%m/%d/%Y') def _vv(v) : return('%.1f%s' % ( v, lbl, )) ax.format_ydata = _vv sa = [ [ datetime.datetime.fromtimestamp(dv[0]), dv[1], ] for dv in sa ] return(sa) def set_size(plt, size) : if size != None : f = plt.gcf() f.set_figwidth( size[0] / float(f.get_dpi())) # turns out, on spring the dpi is 80, but the full image is made using 100 - or something is going on f.set_figheight(size[1] / float(f.get_dpi())) f.set_figwidth( size[0] / 100.0) f.set_figheight(size[1] / 100.0) pass help_str = """ %s (options) gpx_kmz_kml_file(s) Plot point to (+%u) point speed as a function of when the points were and how steep the altitude change is between the points. Note: This program has been tested against "hikified" and "bikified" GPS track files. That is, the points in the tracks have been stripped of redundant points (points in a line between the point's two neighboring points). And, the slope of ascent or descent is determined by comparing the altitude of a point with another point N "steps" ahead in the track. Such N-step logic is sensitive to how the points are distributed in time and space. This program is tuned, therefore, to "hikified/bikified" tracks of someone hiking or biking. Changing the N "step" value with the --step option changes the values that are plotted, as does the --maximize option. To speed up inputting many, many GPS track files, this program can write out a .pickle file with the information and then read the .pickle file in as input instead of reading the original .gpx/.kml/.kmz/etc file(s). The --step and --maximize options have no effect when reading input from a .pickle file. Options: --dds Show or make an image of average/median absolute degrees of slope per day. --dsp Show or make an image of average/median speed per day. --ssl Show or make an image of speed per degrees of slope. --tight Tighten the image margins. --size width hite Set the image size in pixels (does not affect PDF or SVG image files). --label text Set text of the first "word" of the graph titles. --output file_name Put the image(s) to the given file name. If more than one image, add _dds/_dsp/etc. to the base file name. Unless the name already contains a %%t, which is replaced by _dds/_dsp/etc. in any case. The file name extension can be any of the usual image file types. E.g. .png .jpg .bmp .gif .svg etc. Or it can be a .pdf file. If so, and more than one of --dds --dsp --ssl are specified, the images will be put on multiple pages of one .pdf file. --output file_name.json Or, if the --output file name has a .json extension, write the data from the input files to that file for later input instead of the GPS track files. --output file_name.pickle_file Or, if the --output file name has a .pickle extension, write the data from the input files to that file for later input instead of the GPS track files. --median Show median values rather than average. --show Show the plots on-screen. --json Print JSON output. --step N Set how many points ahead the 'next' point is. (default: %u) --maximize Combine the furthest away and furthest up/down in the --step points as the 'next' point. (Warning: Bogus data alert!) --exclude file_name Exclude any matching file name from input. Multiple --excludes allowed. --verbose Increase the verbosity level. """ # # # Main program # # if __name__ == '__main__' : import TZCommandLineAtFile program_name = sys.argv.pop(0) TZCommandLineAtFile.expand_at_sign_command_line_files(sys.argv) dds = False dsp = False ssl = False tight = False label = "" ofile_name = None show_average = True show_plot = False show_json = False json_name = None pickle_name = None size = None STEP = 10 step = None maximize = False excludes = [] verbose = 0 while True : oi = tzlib.array_find(sys.argv, [ "--help", "-h", "-?", "/h", "/?", ] ) if oi < 0 : break del sys.argv[oi] print >> sys.stderr, help_str % ( program_name, STEP, STEP, ) sys.exit(254) while True : oi = tzlib.find_argi(sys.argv, [ "--verbose", "-v", ] ) if oi < 0 : break del sys.argv[oi] verbose += 1 while True : oi = tzlib.find_argi(sys.argv, [ "--dds", ] ) if oi < 0 : break del sys.argv[oi] dds = True while True : oi = tzlib.find_argi(sys.argv, [ "--dsp", ] ) if oi < 0 : break del sys.argv[oi] dsp = True while True : oi = tzlib.find_argi(sys.argv, [ "--ssl", ] ) if oi < 0 : break del sys.argv[oi] ssl = True while True : oi = tzlib.find_argi(sys.argv, [ "--tight", ] ) if oi < 0 : break del sys.argv[oi] tight = True while True : oi = tzlib.find_argi(sys.argv, [ "--size", ] ) if oi < 0 : break del sys.argv[oi] w = max(10, int(sys.argv.pop(oi))) h = max(10, int(sys.argv.pop(oi))) size = [ w, h ] while True : oi = tzlib.find_argi(sys.argv, [ "--label", ] ) if oi < 0 : break del sys.argv[oi] label = sys.argv.pop(oi) while True : oi = tzlib.find_argi(sys.argv, [ "--output", "--out", "--out_file", "--out-file", "--outfile", "--ofile_name", "--ofile-name", "--ofilename", "-o", ] ) if oi < 0 : break del sys.argv[oi] ofile_name = sys.argv.pop(oi) while True : oi = tzlib.find_argi(sys.argv, [ "--step", ] ) if oi < 0 : break del sys.argv[oi] step = max(1, int(sys.argv.pop(oi))) while True : oi = tzlib.find_argi(sys.argv, [ "--maximize", ] ) if oi < 0 : break del sys.argv[oi] maximize = True while True : oi = tzlib.find_argi(sys.argv, [ "--median", "-m", ] ) if oi < 0 : break del sys.argv[oi] show_average = False while True : oi = tzlib.find_argi(sys.argv, [ "--plot", "--show", ] ) if oi < 0 : break del sys.argv[oi] show_plot = True while True : oi = tzlib.find_argi(sys.argv, [ "--exclude", "--excludes", ] ) if oi < 0 : break del sys.argv[oi] excludes.append(sys.argv.pop(oi)) while True : oi = tzlib.find_argi(sys.argv, [ "--json", "-j", ] ) if oi < 0 : break del sys.argv[oi] show_json = True if ofile_name : if os.path.splitext(ofile_name)[1].lower() == '.json' : json_name = ofile_name ofile_name = None elif os.path.splitext(ofile_name)[1].lower() == '.pickle' : pickle_name = ofile_name ofile_name = None elif (dds + dsp + ssl > 1) and (ofile_name.find('%t') < 0) : fn, ext = os.path.splitext(ofile_name) ofile_name = fn + '%t' + ext pass tracks = [] speeds = None # when not None, this is a dictionary keyed by track file name, valued by an array of arrays with [ angle, speed, when, track file_name ] try : while len(sys.argv) : fn = sys.argv.pop(0) if os.path.splitext(fn)[1].lower() == '.json' : speeds = json.loads(tzlib.read_whole_text_file(fn)) if len(tracks) : print "Tracks already read and you are trying to read their data from %s!" % fn sys.exit(108) if len(sys.argv) : print "More command line things after the .json input file: %s!" % " ".join(sys.argv) sys.exit(109) if step != None : print "Warning: --step has no effect. It has already been used to make the contents of %s" % fn if maximize : print "Warning: --maximize has no effect. It has already been used to make the contents of %s" % fn if len(excludes) : print "Excludes don't affect data read from %s" % fn break elif os.path.splitext(fn)[1].lower() == '.pickle' : speeds = tzlib.unpickle_file(fn) if len(tracks) : print "Tracks already read and you are trying to read their data from %s!" % fn sys.exit(108) if len(sys.argv) : print "More command line things after the .pickle input file: %s!" % " ".join(sys.argv) sys.exit(109) if step != None : print "Warning: --step has no effect. It has already been used to make the contents of %s" % fn if maximize : print "Warning: --maximize has no effect. It has already been used to make the contents of %s" % fn if len(excludes) : print "Excludes don't affect data read from %s" % fn break fns = tzlib.ambiguous_file_list(fn, do_sub_dirs = False) for fnn in fns : mtc = None for fnnn in excludes : if fnmatch.fnmatch(fnn, fnnn) or ((not os.path.dirname(fnnn)) and fnmatch.fnmatch(os.path.basename(fn), fnnn)) : mtc = fnnn break pass if not mtc : if verbose > 1 : print fnn ta = tz_gps.tracks_from_file(fnn) # note: the file name can be converted to a file name in metadata in a .gpx file - an abspath file name, in practice at least one time: muriel_fresh_and_easy_01x.gpx if not ta : print "%sNo tracks in %s!" % ( ((verbose > 1) and " ") or "", fnn, ) else : tracks += ta pass elif verbose > 2 : print "Exclude: %s by %s" % ( fn, mtc, ) pass pass pass except KeyboardInterrupt : pass if speeds is None : if not len(tracks) : print "No tracks found!" sys.exit(101) start_when = min([ t.points[ 0].when for t in tracks ]) end_when = max([ t.points[-1].when for t in tracks ]) if step is None : step = STEP speeds = {} for t in tracks : fn = t.file_name if fn not in speeds : speeds[fn] = [] for pi, p in enumerate(t.points[0:-step]) : if maximize : mxa = 0 xxa = 0.0 mxd = 0 mxp = p for pii in xrange(pi + 1, pi + step + 1) : pp = t.points[pii] d = abs(pp.altitude - p.altitude) if mxa < d : mxa = d xxa = pp.altitude d = p.flat_distance_from(pp) if mxd < d : mxd = d mxp = pp pass pp = copy.deepcopy(mxp) pp.altitude = xxa else : pp = t.points[pi + step] a = math.degrees(math.atan2(pp.altitude - p.altitude, p.flat_distance_from(pp) * latlon.metersPerNauticalMile)) speeds[fn].append([ a, p.geo_speed_from(pp), p.when, fn, ]) if not len(speeds[fn]) : del(speeds[fn]) pass pass else : start_when = min([ min([ sss[2] for sss in ss ]) for ss in speeds.values() ]) end_when = max([ max([ sss[2] for sss in ss ]) for ss in speeds.values() ]) if verbose : print "%u tracks, start %s to end %s for %u days, %u points" % ( len(speeds), time.asctime(time.localtime(start_when)), time.asctime(time.localtime(end_when)), 1 + ((end_when - start_when) / A_DAY), point_count(speeds), ) sys.argv.insert(0, program_name) # make matplotlib's logic not mess up if json_name : tzlib.write_whole_text_file(json_name, json.dumps(speeds, indent = 4)) elif pickle_name : if not tzlib.pickle_file(pickle_name, speeds, proto = sys.maxsize) : print "Failed to write %s!" % pickle_name sys.exit(105) pass elif verbose and (not dds) and (not dsp) and (not ssl) and show_json : # sometime do the original 3d contour map interpolated with RBF interpolation, griddata, or meshgrid or whatever for a in speeds : a.sort(lambda a, b : cmp(a[0], b[0]) or cmp(a[2], b[2])) # sort by angle, per when if show_json : print json.dumps(speeds, indent = 4) pass pass if ofile_name and (dds + dsp + ssl > 1) and (os.path.splitext(ofile_name)[1].lower() == '.pdf') : ofile_name = ofile_name.replace('%t', '') pdf = PdfPages(ofile_name) else : pdf = None if dds : if len(speeds) < 2 : print "Cannot --dds graph a single track." else : sa = [] for a in speeds.values() : t = a[0][2] if show_average : sa.append([ t, sum([ abs(asp[0]) for asp in a ]) / len(a), ]) g_typ = "Average" else : sa.append([ t, tzlib.median([ abs(asp[0]) for asp in a ]), ]) g_typ = "Median" if (verbose > 3) and (sa[-1][1] > 20) : print "High angle track: Middle point's abs angle, speed, when, file name:", a[len(a) / 2] pass if show_json : print json.dumps(sa, indent = 4) if show_plot or ofile_name : fig = plt.figure() sa = plot_years(plt, sa, start_when, end_when, u'\u00b0') plt.plot([ a[0] for a in sa ], [ a[1] for a in sa ], 'ro') if plt.ylim()[0] > 0 : plt.ylim(0.0, plt.ylim()[1]) ttl = "%sPer-Track %s Absolute Degree of Slope" % ( label, g_typ, ) plt.title(ttl) plt.gcf().canvas.set_window_title(ttl) plt.xlabel("Date") plt.ylabel("Absolute Degree of Slope") plt.grid(True) if tight : plt.tight_layout() set_size(plt, size) if ofile_name : if pdf : pdf.savefig(fig) else : plt.savefig(ofile_name.replace('%t', '_dds')) pass if show_plot : plt.show() pass pass pass if dsp : if len(speeds) < 2 : print "Cannot --dsp graph a single track." else : sa = [] for a in speeds.values() : t = a[0][2] if show_average : if False : a = list(a) a.sort(lambda a, b : cmp(abs(a[0]), abs(b[0])) or cmp(a[1], b[1]) or cmp(a[2], b[2])) a = a[:10] sa.append([ t, sum([ asp[1] for asp in a ]) / len(a), ]) g_typ = "Average" else : sa.append([ t, tzlib.median([ asp[1] for asp in a ]), ]) g_typ = "Median" pass if show_json : print json.dumps(sa, indent = 4) if show_plot or ofile_name : fig = plt.figure() sa = plot_years(plt, sa, start_when, end_when, ' kph') plt.plot([ aa[0] for aa in sa ], [ aa[1] for aa in sa ], 'ro') if plt.ylim()[0] > 0 : plt.ylim(0.0, plt.ylim()[1]) ttl = "%sPer-Track %s KPH Speeds" % ( label, g_typ, ) plt.title(ttl) plt.gcf().canvas.set_window_title(ttl) plt.xlabel("Date") plt.ylabel("KPH Speed") plt.grid(True) # print "@@@@", plt.gcf().get_size_inches(), plt.gcf().get_dpi() if tight : plt.tight_layout() set_size(plt, size) if ofile_name : if pdf : pdf.savefig(fig) else : plt.savefig(ofile_name.replace('%t', '_dsp')) pass if show_plot : plt.show() pass pass pass if ssl : sa = [] for a in speeds.values() : if a != None : sa += a pass sa.sort(lambda a, b : cmp(a[0], b[0]) or cmp(a[1], b[1]) or cmp(a[2], b[2])) if show_json : print json.dumps(sa, indent = 4) if show_plot or ofile_name : fig = plt.figure() def _xv(v) : return(u'%.1f\u00b0' % v) plt.gca().format_xdata = _xv def _yv(v) : return(u'%.1f kph' % v) plt.gca().format_ydata = _yv plt.plot([ aa[0] for aa in sa ], [ aa[1] for aa in sa ], 'ro') if plt.ylim()[0] > 0 : plt.ylim(0.0, plt.ylim()[1]) ttl = "%sKPH Speeds per Slope Degrees" % label plt.title(ttl) plt.gcf().canvas.set_window_title(ttl) # of plt.gca().get_title() plt.xlabel("Slope in Degrees") plt.ylabel("KPH Speed") plt.grid(True) if tight : plt.tight_layout() set_size(plt, size) if ofile_name : if pdf : pdf.savefig(fig) else : plt.savefig(ofile_name.replace('%t', '_ssl')) pass if show_plot : plt.show() pass pass if pdf : d = pdf.infodict() d['Title'] = os.path.basename(ofile_name) un = tzlib.get_full_user_name() if un : d['Author'] = un d['Subject'] = "Graphs of GPS Tracks' Speed and Slope" d['Keywords'] = "GPS,speed,slope,angle" d['CreationDate'] = datetime.datetime.utcfromtimestamp(end_when) d['ModDate'] = datetime.datetime.today() d['Producer'] = os.path.splitext(os.path.basename(__file__))[0] + '.py' pdf.close() pass # # # eof