package net.sourceforge.gpstools;
/* gpsdings
* Copyright (C) 2006 Moritz Ringler
* $Id: TrackAnalyzer.java 445 2011-01-17 21:26:27Z ringler $
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
import java.io.BufferedInputStream;
import java.io.BufferedOutputStream;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.InputStream;
import java.io.OutputStream;
import java.io.OutputStreamWriter;
import java.io.PrintStream;
import java.text.DateFormat;
import java.text.SimpleDateFormat;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Date;
import java.util.Enumeration;
import java.util.List;
import java.util.SortedSet;
import java.util.Properties;
import net.sourceforge.gpstools.gpx.*;
import net.sourceforge.gpstools.kml.GpxTrackWriter;
import net.sourceforge.gpstools.kml.QuantityTrackWriter;
import net.sourceforge.gpstools.plot.SvgPlot;
import net.sourceforge.gpstools.utils.*;
import net.sourceforge.gpstools.math.*;
import org.apache.commons.cli.CommandLine;
import org.apache.commons.cli.Options;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
import org.apache.commons.math.analysis.PolynomialSplineFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.UnivariateRealSolver;
import org.apache.commons.math.analysis.UnivariateRealSolverFactory;
import org.apache.commons.math.analysis.SplineInterpolator;
import org.apache.commons.math.complex.Complex;
import org.xml.sax.SAXException;
public class TrackAnalyzer extends GPSDings {
private final Interpolator interpolator = new Interpolator();
private float speedCutoff = 1.0f;
private boolean stats = true;
private boolean usUnits = false;
private final static float SPEED_UNIT = 3600000.0f;
private WptType trackOrigin;
private GpxType mygpx = new Gpx();
public enum Quantity {
TimeTravelled("Time Travelled (min)"), DistanceTravelled(
"Distance Travelled (km)"), Latitude("Latitude (deg)"), Longitude(
"Longitude (deg)"), Elevation("Elevation (m)"), HorizontalSpeed(
"Horizontal Speed (km/h)"), VerticalSpeed(
"Vertical Speed (m/h)"), Slope("Slope (%)"), Bearing(
"Bearing (deg)");
private final String longname;
Quantity(String longname) {
this.longname = longname;
}
public String longName() {
return longname;
}
}
private double totalAscent = 0.0;
private double totalDescent = 0.0;
private UnivariateRealSolverFactory sf;
protected SplineInterpolator spline = new SplineInterpolator();
public TrackAnalyzer() {
// sole constructor
}
public int addAllTracks(GpxType gpx) {
Enumeration<? extends Trk> trks = gpx.enumerateTrk();
while (trks.hasMoreElements()) {
mygpx.addTrk(trks.nextElement());
}
return interpolator.addAllTracks(gpx);
}
public void addTrack(TrkType trk) {
interpolator.addTrack(trk);
}
public void addTrackSegment(TrksegType seg) {
interpolator.addTrackSegment(seg);
}
/** decimal degrees */
public double getLat(Date d) throws FunctionEvaluationException {
return interpolator.fLati.value(d.getTime());
}
/** decimal degrees */
public double getLon(Date d) throws FunctionEvaluationException {
return interpolator.fLongi.value(d.getTime());
}
/** m */
public double getEle(Date d) throws FunctionEvaluationException {
return (interpolator.fEle == null) ? 0.0 : interpolator.fEle.value(d
.getTime());
}
/** km */
public double getDist(Date d) throws FunctionEvaluationException {
long dd = d.getTime();
return interpolator.fDist.value(dd);
}
/** km/h */
public double getHorizontalSpeed(Date d) throws FunctionEvaluationException {
return Math
.max(0, interpolator.fHSpeed.value(d.getTime()) * SPEED_UNIT);
}
/** m/h */
public double getVerticalSpeed(Date d) throws FunctionEvaluationException {
return (interpolator.fVSpeed == null) ? 0.0 : interpolator.fVSpeed
.value(d.getTime()) * 3600000;
}
/** m/km */
public double getSlope(Date d) throws FunctionEvaluationException {
double dtime = d.getTime();
return interpolator.fVSpeed.value(dtime)
/ interpolator.fHSpeed.value(dtime) / 10;
}
/** deg **/
public double getBearing(Date d) throws FunctionEvaluationException {
return interpolator.fBearing.value(d.getTime());
}
public double getQuantity(Date d, Quantity q)
throws FunctionEvaluationException {
switch (q) {
case TimeTravelled:
return (d.getTime() - getOrigin().getTime().getTime()) / 60000.;
case DistanceTravelled:
return getDist(d);
case Latitude:
return getLat(d);
case Longitude:
return getLon(d);
case Elevation:
return getEle(d);
case HorizontalSpeed:
return getHorizontalSpeed(d);
case VerticalSpeed:
return getVerticalSpeed(d);
case Slope:
return getSlope(d);
case Bearing:
return getBearing(d);
default:
throw new Error("Unknown quantity " + q);
}
}
private double unitConversion(Quantity q, double val) {
double result;
if (usUnits) {
switch (q) {
case DistanceTravelled: // falls through
case HorizontalSpeed:
result = kmToMi(val);
break;
case Elevation: // falls through
case VerticalSpeed:
result = mToFt(val);
break;
default:
result = val;
}
} else {
result = val;
}
return result;
}
private String getAxisTitle(Quantity q) {
String result;
if (usUnits) {
switch (q) {
case DistanceTravelled:
result = "Distance Travelled (mi)";
break;
case HorizontalSpeed:
result = "Horizontal Speed (mi/h)";
break;
case Elevation:
result = "Elevation (ft)";
break;
case VerticalSpeed:
result = "Vertical Speed (ft/h)";
break;
default:
result = q.longName();
}
} else {
result = q.longName();
}
return result;
}
public void plot(Quantity x, Quantity y, OutputStream os, String title,
Wpt[] pts) throws IOException {
OutputStreamWriter ow = new OutputStreamWriter(os, "utf-8");
try {
SvgPlot plot = getPlot(x, y, title, pts);
configurePlot(plot, x, y);
plot.plot(ow);
ow.flush();
} catch (FunctionEvaluationException ex) {
// should never happen.
throw new Error(ex);
}
}
public SvgPlot getPlot(Quantity x, Quantity y, String title, Wpt[] pts)
throws FunctionEvaluationException {
SvgPlot plot = new SvgPlot();
plot.setTitle(SvgPlot.PlotDesignation.GLOBAL, title);
for (double[] xx : getSupportPoints()) {
final int n = xx.length;
double[] xdata = new double[n];
double[] ydata = new double[n];
for (int i = 0; i < n; i++) {
Date d = new Date((long) xx[i]);
xdata[i] = unitConversion(x, getQuantity(d, x));
ydata[i] = unitConversion(y, getQuantity(d, y));
}
plot.addData(xdata, ydata);
}
for (Wpt wpt : pts) {
for (Date d : timeWpt(wpt, 0.01)) {
try {
plot.addPoint(unitConversion(x, getQuantity(d, x)),
unitConversion(y, getQuantity(d, y)), wpt.getName());
} catch (Exception err) {
System.err.println("Warning: Cannot plot Wpt:");
System.err.println(wpt.getName());
System.err.println(wpt.getLon());
System.err.println(wpt.getLat());
DateFormat df = new SimpleDateFormat(GPSDings.ISO_DATE
+ "'Z'");
df.setTimeZone(GPSDings.UTC);
System.err.println(df.format(d));
System.err.println(err);
}
}
}
plot.setTitle(SvgPlot.PlotDesignation.X, getAxisTitle(x));
plot.setTitle(SvgPlot.PlotDesignation.Y, getAxisTitle(y));
return plot;
}
private void configurePlot(SvgPlot plot, Quantity x, Quantity y) {
String pFile = System.getProperty(getClass().getName() + ".properties");
if (pFile != null) {
InputStream in = null;
try {
Properties p = new Properties();
in = new FileInputStream(pFile);
in = new BufferedInputStream(in);
p.load(in);
String prop;
prop = p.getProperty("Scale");
if (prop != null) {
plot.setScale(Float.parseFloat(prop));
}
prop = p.getProperty("y0");
int y0 = (prop == null) ? -1 : Integer.parseInt(prop);
prop = p.getProperty("x0");
int x0 = (prop == null) ? -1 : Integer.parseInt(prop);
if (x0 > 0 && y0 > 0) {
plot.setOrigin(x0, y0);
}
prop = p.getProperty(x.name());
if (prop != null) {
plot.setTitle(SvgPlot.PlotDesignation.X, prop);
}
prop = p.getProperty(y.name());
if (prop != null) {
plot.setTitle(SvgPlot.PlotDesignation.Y, prop);
}
} catch (Exception ex) {
message(ex);
message("Warning: cannot load SvgPlot configuration from "
+ pFile);
} finally {
if (in != null) {
try {
in.close();
} catch (IOException ex) {
System.err.println(ex);
}
}
}
}
}
/**
* @throws IllegalArgumentException
* Argument must lie between 0 and 20
**/
public void setSmoothing(int i) {
if (i < 0 || i > 20) {
throw new IllegalArgumentException(
"Argument must lie between 0 and 20");
}
interpolator.smooth = i;
}
public void setSpeedCutoff(float f) {
speedCutoff = f;
}
public void setPrintStats(boolean b) {
stats = b;
}
public int getSmoothing() {
return interpolator.smooth;
}
private synchronized void maybeSetOrigin(WptType candidate) {
if (trackOrigin == null
|| TrackInterpolator.compareByDateTime(candidate, trackOrigin) < 0) {
trackOrigin = candidate;
}
}
public synchronized WptType getOrigin() {
return trackOrigin;
}
public Date[] timeWpt(WptType p, double maxDist)
throws FunctionEvaluationException {
/* check if the waypoint already has a time */
if (p.getTime() != null) {
return new Date[] { p.getTime() };
}
List<Date> result = new ArrayList<Date>();
double pLat = p.getLat().doubleValue();
double pLon = p.getLon().doubleValue();
/* check which of the track points are closest */
double[][] dd = getSupportPoints();
for (double[] d : dd) {
double dist;
for (int i = 0, n = d.length; i < n; i++) {
/* calculate distance of trackpoint i */
dist = GpsMath.kmDistance(interpolator.fLati.value(d[i]),
interpolator.fLongi.value(d[i]), pLat, pLon);
/* if trackpoint i is closer than our threshold distance */
if (dist <= maxDist) {
int iclosest = i;
double dclosest = dist;
/*
* check whether the following trackpoints are even closer.
*/
for (i++; i < n; i++) {
dist = GpsMath.kmDistance(
interpolator.fLati.value(d[i]),
interpolator.fLongi.value(d[i]), pLat, pLon);
if (dist < dclosest) {
iclosest = i;
dclosest = dist;
} else if (dist > maxDist) {
break;
}
}
result.add(new Date((long) d[iclosest]));
}
}
}
return result.toArray(new Date[result.size()]);
}
public double[][] getSupportPoints() {
PieceWiseFunction.ClosedInterval[] domain = interpolator.fDist
.getDomain();
final int n = domain.length;
double[][] result = new double[n][];
try {
for (int i = 0; i < n; i++) {
result[i] = ((LinearInterpolator.LinearInterpolation) interpolator.fDist
.getSegmentFor(domain[i].xmin)).getKnots();
}
} catch (RuntimeException ex) {
throw ex;
} catch (Exception ex) {
throw new RuntimeException(ex);
}
return result;
}
public void printKML(Quantity q, OutputStream out, String name)
throws FunctionEvaluationException, SAXException, IOException {
GpxTrackWriter qtw = QuantityTrackWriter
.getInstance(out, name, q, this);
qtw.writeGpx(mygpx);
qtw.endDocument();
}
/**
* Generates output of the specified type.
*
* @param out
* the type of output to produce.
* @throws IOException
* if an I/O error occurs.
*/
@Override
public void print(Output out) throws IOException {
if (out == Output.Text) {
File f = output.get(out);
PrintStream os = (f == null) ? System.out : new PrintStream(
new BufferedOutputStream(new FileOutputStream(f)));
printText(os);
os.flush();
os.close();
} else {
throw new IllegalArgumentException("Output mode " + out
+ "is not supported by TrackAnalyzer.");
}
}
@Override
public void print() throws IOException {
super.print();
if (stats) {
printStats(System.err);
}
}
private void printStats(PrintStream ps) throws IOException {
try {
long timeNotMoving = 0l;
final double[][] supp = getSupportPoints();
final int nseg = supp.length;
double lastX = Double.NaN;
double lastT = Double.NaN;
for (int i = 0; i < nseg; i++) {
double[] seg = supp[i];
final int seglength = seg.length;
for (int j = 0; j < seglength; j++) {
double time = seg[j];
double dist = getDist(new Date((long) time));
// see if we have moved
if (!Double.isNaN(lastX)) {
double deltaX = dist - lastX;
double deltaT = time - lastT;
if (deltaX / deltaT * SPEED_UNIT < speedCutoff) {
timeNotMoving += (long) deltaT;
}
}
lastX = dist;
lastT = time;
}
}
final long last = (long) lastT;
final long totalTime = last - getOrigin().getTime().getTime();
final long timeMoving = totalTime - timeNotMoving;
double distKm = lastX;
long time = totalTime;
ps.printf("TimeTravelled\t%02d:%02d:%02.0f\n", (time / 3600000),
(time / 60000) % 60, (time / 1000.) - (time / 60000 * 60));
time = timeNotMoving;
ps.printf("TimeNotMoving\t%02d:%02d:%02.0f\n", (time / 3600000),
(time / 60000) % 60, (time / 1000.) - (time / 60000 * 60));
time = timeMoving;
ps.printf("TimeMoving\t%02d:%02d:%02.0f\n", (time / 3600000),
(time / 60000) % 60, (time / 1000.) - (time / 60000 * 60));
String km = "km";
String m = "m";
double totalAsc = totalAscent;
double totalDesc = totalDescent;
if (usUnits) {
distKm = kmToMi(distKm);
totalAsc = mToFt(totalAsc);
totalDesc = mToFt(totalDesc);
km = "mi";
m = "ft";
}
ps.printf("DistanceTravelled\t%5.3f %s\n", distKm, km);
ps.printf("TotalAscent\t%5.3f %s\n", totalAsc, m);
ps.printf("TotalDescent\t%5.3f %s\n", totalDesc, m);
time = totalTime;
ps.printf("AverageSpeed\t%5.2f %s/h\n", distKm / (time / 3600000.),
km);
time = timeMoving;
ps.printf("AverageSpeedWhileMoving\t%5.2f %s/h\n", distKm
/ (time / 3600000.), km);
ps.printf("(Speeds less than %2.1f km/h were treated as zero)\n",
speedCutoff);
} catch (FunctionEvaluationException ex) {
// should never happen.
throw new Error(ex);
}
}
public void setUseUSUnits(boolean b) {
usUnits = b;
}
public boolean getUseUSUnits() {
return usUnits;
}
private double mToFt(double mval) {
return 3.2808399 * mval;
}
private double kmToMi(double kmval) {
return 0.62137119 * kmval;
}
private void printText(PrintStream ps) throws IOException {
if (getOrigin() == null) {
return;
}
double timezero = getOrigin().getTime().getTime();
if (usUnits) {
ps.println("Date Time/min Lat Lon Dist/mi Horizontal Speed/(mi/h) Bearing/deg Elevation/ft Vertical Speed/(ft/h) Slope/%"
.replaceAll(" ", "\t"));
} else {
ps.println("Date Time/min Lat Lon Dist/km Horizontal Speed/(km/h) Bearing/deg Elevation/m Vertical Speed/(m/h) Slope/%"
.replaceAll(" ", "\t"));
}
GpsFormat format = GpsFormat.getInstance();
DateFormat df = new SimpleDateFormat(GPSDings.ISO_DATE + "'Z'");
df.setTimeZone(GPSDings.UTC);
try {
for (double[] xx : getSupportPoints()) {
ps.println("#new segment");
for (double x : xx) {
Date d = new Date((long) x);
ps.print(df.format(d));
ps.printf("\t%5.2f", (x - timezero) / 60000.);
ps.print("\t" + format.asLatLon(getLat(d)));
ps.print("\t" + format.asLatLon(getLon(d)));
ps.printf("\t%5.3f", usUnits ? kmToMi(getDist(d))
: getDist(d));
ps.printf("\t%3.1f",
usUnits ? kmToMi(getHorizontalSpeed(d))
: getHorizontalSpeed(d));
ps.printf("\t%3.1f", getBearing(d));
ps.printf("\t%4.1f", usUnits ? mToFt(getEle(d)) : getEle(d));
ps.printf("\t%3.1f", usUnits ? mToFt(getVerticalSpeed(d))
: getVerticalSpeed(d));
ps.printf("\t%3.1f", getSlope(d));
ps.println();
}
}
} catch (FunctionEvaluationException ex) {
// should never happen.
throw new Error(ex);
}
}
public static void main(String[] argv) throws Exception {
(new AnalyzerCommandLine(argv)).execute();
}
private static Complex[] applyHanningWindow(Complex[] z, int windowWidth) {
int width = windowWidth;
final int n = z.length;
if (width < 2) {
if (n >= 2) {
width = 2;
} else {
return z;
}
}
Complex[] filtered = new Complex[n];
int j = 0;
for (final int limit = width / 2; j < limit; j++) {
filtered[j] = z[j].multiply(new Complex(0.53836 + 0.46164 * Math
.cos(2 * Math.PI * j / width), 0.0));
}
for (final int limit = n - width / 2; j < limit; j++) {
filtered[j] = Complex.ZERO;
}
for (int i = -width / 2; j < n; j++) {
filtered[j] = z[j].multiply(new Complex(0.53836 + 0.46164 * Math
.cos(2 * Math.PI * i / width), 0.0));
i++;
}
return filtered;
}
private class Interpolator extends TrackInterpolator {
final PieceWiseFunction fDist = new PieceWiseFunction() {
@Override
public double value(double x) throws FunctionEvaluationException {
double offset = 0.0;
SegmentFunction last = null;
for (SegmentFunction f : elements) {
if (last != null) {
offset += GpsMath.kmDistance(getWpt(new Date(
(long) last.domain.xmax)), getWpt(new Date(
(long) f.domain.xmin)));
}
if (f.domain.contains(x)) {
return f.func.value(x) + offset;
}
offset += f.func.value(f.domain.xmax);
last = f;
}
throw new FunctionEvaluationException(x,
"Argument outside domain: " + x);
}
};
final PieceWiseFunction fHSpeed = new PieceWiseFunction();
final PieceWiseFunction fVSpeed = new PieceWiseFunction();
final PieceWiseFunction fBearing = new PieceWiseFunction();
volatile int smooth = 1;
public Interpolator() {
// no extra initialization.
}
@Override
protected synchronized void interpolateSegment(SortedSet<WptType> points) {
final int n = points.size();
if (n == 0) {
return;
}
super.interpolateSegment(points);
double[] time = new double[n];
double[] dist = new double[n];
final WptType firstWpt = points.iterator().next();
WptType lastWpt = firstWpt;
double length = 0.0;
int i = 0;
points.size();
for (WptType wpt : points) {
time[i] = wpt.getTime().getTime();
length += GpsMath.kmDistance(lastWpt, wpt);
lastWpt = wpt;
dist[i] = length;
i++;
}
lastWpt = null;
final double tmin = time[0];
final double tmax = time[n - 1];
try {
fDist.addSegmentFunction(tmin, tmax,
super.interpolator.interpolate(time, dist));
} catch (Exception wonthappen) {
throw new Error(wonthappen);
}
maybeSetOrigin(firstWpt);
try {
fHSpeed.addSegmentFunction(tmin, tmax, smooth(time, dist)
.derivative());
fVSpeed.addSegmentFunction(tmin, tmax,
smooth(fEle.getSegmentFor(tmin), tmin, tmax, n)
.derivative());
final DifferentiableUnivariateRealFunction latseg = smooth(
fLati.getSegmentFor(tmin), tmin, tmax, n);
fBearing.addSegmentFunction(tmin, tmax, new BearingFunction(
smooth(fLongi.getSegmentFor(tmin), tmin, tmax, n)
.derivative(), latseg.derivative(), latseg));
/* Calculate total ascent and total descent */
List<Double> vpoints = getZeros(fVSpeed, time);
vpoints.add(tmax);
double z0 = fEle.value(tmin);
for (Double d : vpoints) {
double zi = fEle.value(d);
double deltaz = zi - z0;
if (deltaz > 0) {
totalAscent += deltaz;
} else {
totalDescent -= deltaz;
}
z0 = zi;
}
} catch (Exception ex) {
ex.printStackTrace();
throw new RuntimeException(ex);
}
}
private List<Double> getZeros(UnivariateRealFunction f, double[] x)
throws FunctionEvaluationException, ConvergenceException {
List<Double> zeros = new ArrayList<Double>();
final int n = x.length;
if (x.length != 0) {
synchronized (this) {
if (sf == null) {
sf = UnivariateRealSolverFactory.newInstance();
}
}
UnivariateRealSolver solver = sf.newDefaultSolver(f);
double sign = Math.signum(f.value(x[0]));
for (int i = 1; i < n; i++) {
double isig = Math.signum(f.value(x[i]));
if (isig != sign) {
zeros.add(solver.solve(x[i - 1], x[i]));
sign = isig;
}
}
}
return zeros;
}
private DifferentiableUnivariateRealFunction smooth(double[] x,
double[] y) throws org.apache.commons.math.MathException {
return smooth(spline.interpolate(x, y), x[0], x[x.length - 1],
x.length);
}
private DifferentiableUnivariateRealFunction smooth(
UnivariateRealFunction f, double xmin, double xmax,
int minNumPoints) throws org.apache.commons.math.MathException {
if (getSmoothing() == 0
&& f instanceof DifferentiableUnivariateRealFunction) {
return (DifferentiableUnivariateRealFunction) f;
}
/*
* Sample the function at samplenum equidistant points, where
* samplenum is the smallest power of two that is greater than the
* minNumPoints.
*/
int samplenum = 2;
int len = minNumPoints;
while ((len /= 2) != 0) {
samplenum *= 2;
}
final double[] samples = FastFourierTransformer.sample(f, xmin,
xmax, samplenum);
if (getSmoothing() == 0) {
double[] x = linspace(xmin, xmax, samplenum);
return (PolynomialSplineFunction) spline
.interpolate(x, samples);
}
// if we have negative numbers we just shift everything
// so we just have positive numbers and don't lose sign
// in the FFT cycle
// Just a temporary fix until I have time to think this through
double min = Double.POSITIVE_INFINITY;
for (int i = 0; i < samplenum; i++) {
if (samples[i] < min) {
min = samples[i];
}
}
if (min < 0) {
for (int i = 0; i < samplenum; i++) {
samples[i] -= min;
}
}
/*
* To avoid oscillations at the start and at the end of the filtered
* data pad the samples array as follows data = { (samples[0])n/2 ,
* samples, (samples[n-1])n/2 } where n = samplenum
*/
final int datanum = samplenum * 2;
final double[] data = new double[datanum];
System.arraycopy(samples, 0, data, samplenum / 2, samplenum);
Arrays.fill(data, 0, samplenum / 2, samples[0]);
Arrays.fill(data, samplenum * 3 / 2, datanum,
samples[samplenum - 1]);
/* Fourier transform the data array. */
final FastFourierTransformer transformer = new FastFourierTransformer();
final Complex[] fouriertransform = transformer.transform(data);
/*
* Apply a Hanning window to the fourier transform. This has the
* effect of a long-pass filter. The width of the window is
* Math.max( 2, datanum/2^getSmooting() );
*/
final int width = datanum / (1 << getSmoothing());
final Complex[] filtered = applyHanningWindow(fouriertransform,
width);
/*
* Inverse fourier transform the filtered dataset and interpolate it
* to obtain the smoothed function.
*/
final Complex[] smoothed = transformer.inversetransform(filtered);
final double[] newy = new double[samplenum];
final double[] newx = new double[samplenum];
final double xstep = (xmax - xmin) / (samplenum - 1);
min = Math.min(min, 0);
for (int i = 0, k = samplenum / 2; i < samplenum; i++) {
newx[i] = xmin + i * xstep;
newy[i] = smoothed[k++].abs() + min;
}
return (PolynomialSplineFunction) spline.interpolate(newx, newy);
}
}
private static class AnalyzerCommandLine extends
AbstractCommandLine<TrackAnalyzer> {
public static final String OPT_SMOOTH = "smoothing";
public static final String OPT_PLOT = "plot";
public static final String OPT_KML = "kml";
public static final String OPT_TEXT = "text";
public static final String OPT_MIN_SPEED = "vmin";
public static final String OPT_NO_STATS = "no-stats";
public static final String OPT_US_UNITS = "feet-and-miles";
private static class PlotParameters {
File output;
Quantity x;
Quantity y;
String title;
}
private PlotParameters[] plots;
private PlotParameters[] kmls;
public AnalyzerCommandLine(String[] argv) {
super(new TrackAnalyzer(), argv);
}
@Override
public void execute() {
try {
// java.util.Locale.setDefault(java.util.Locale.US);
processArguments();
String[] gpxFile = getInput();
if (gpxFile == null || gpxFile.length != 1) {
throw new IllegalArgumentException(
"You must specify exactly one gpx input file.");
}
String fname = "";
InputStream in = System.in;
if (!"-".equals(gpxFile[0])) {
File f = new File(gpxFile[0]);
in = new FileInputStream(f);
fname = f.getName() + ": ";
}
GpxType xgpx = null;
try {
xgpx = GPSDings.readGPX(in);
} finally {
in.close();
}
app.addAllTracks(xgpx);
if (app.getSupportPoints().length == 0) {
app.message("No tracks with time information found. Trying to parameterize by distance.");
Trkpt last = null;
long time = 0l;
for (Trk trk : xgpx.getTrk()) {
for (Trkseg seg : trk.getTrkseg()) {
for (Trkpt pt : seg.getTrkpt()) {
if (last != null) {
time += (long) (GpsMath
.kmDistance(last, pt) * 3600 * 1000);
}
pt.setTime(new Date(time));
last = pt;
}
}
}
app.mygpx = new Gpx();
app.addAllTracks(xgpx);
}
app.print();
if (plots != null) {
for (PlotParameters pp : plots) {
OutputStream os = new FileOutputStream(pp.output);
try {
os = new BufferedOutputStream(os);
if (pp.title == null) {
pp.title = fname + pp.y + " vs. " + pp.x;
}
app.plot(pp.x, pp.y, os, pp.title, xgpx.getWpt());
} finally {
os.close();
}
}
}
if (kmls != null) {
for (PlotParameters pp : kmls) {
OutputStream os = new FileOutputStream(pp.output);
try {
os = new BufferedOutputStream(os);
if (pp.title == null) {
pp.title = fname + pp.x;
}
app.printKML(pp.x, os, pp.title);
} finally {
os.close();
}
}
}
} catch (Exception ex) {
app.handleException(null, ex);
System.exit(1);
}
}
@Override
public void printHelp(PrintStream out) {
try {
super.printHelp(out);
cat("trackanalyzer.txt");
} catch (Exception ex) {
throw new Error(ex);
}
}
@Override
protected Options createOptions() {
Options options = super.createOptions();
options.addOption(makeOption(OPT_TEXT, File.class, 1, "file"));
options.addOption(makeOption(OPT_SMOOTH, Integer.class, 1, "int"));
options.addOption(makeOption(OPT_PLOT, String.class, 1, "plotlist"));
options.addOption(makeOption(OPT_KML, String.class, 1, "paramlist"));
options.addOption(makeOption(OPT_MIN_SPEED, Float.class, 1, "f"));
options.addOption(makeOption(OPT_NO_STATS, Boolean.class, 0, null));
options.addOption(makeOption(OPT_US_UNITS, Boolean.class, 0, null));
return options;
}
@Override
protected CommandLine processOptions(Options options)
throws org.apache.commons.cli.ParseException, IOException,
java.text.ParseException {
CommandLine cl = super.processOptions(options);
char c;
/* smoothing option */
opt = OPT_SMOOTH;
c = opt.charAt(0);
if (cl.hasOption(c)) {
app.setSmoothing(Integer.parseInt(cl.getOptionValue(c)));
}
/* plot option */
opt = OPT_PLOT;
c = opt.charAt(0);
if (cl.hasOption(c)) {
String plotlist = cl.getOptionValue(c);
String[] plotparam = plotlist.split(";", -1);
List<PlotParameters> lpp = new ArrayList<PlotParameters>();
for (String paramlist : plotparam) {
String[] params = paramlist.split("\\|", -1);
if (params.length < 3 || params.length > 4) {
throw new IllegalArgumentException(
"Illegal plot specification "
+ paramlist
+ "\nPlots must be specified as Qx|Qy|File|Title"
+ "\nMultiple plots can be separated by ';'.");
}
PlotParameters pp = new PlotParameters();
pp.x = Enum.valueOf(Quantity.class, params[0]);
pp.y = Enum.valueOf(Quantity.class, params[1]);
pp.output = new File(params[2]);
if (params.length == 4) {
pp.title = params[3];
}
lpp.add(pp);
}
plots = lpp.toArray(new PlotParameters[lpp.size()]);
}
/* kml option */
opt = OPT_KML;
c = opt.charAt(0);
if (cl.hasOption(c)) {
String plotlist = cl.getOptionValue(c);
String[] plotparam = plotlist.split(";", -1);
List<PlotParameters> lpp = new ArrayList<PlotParameters>();
for (String paramlist : plotparam) {
String[] params = paramlist.split("\\|", -1);
if (params.length < 2 || params.length > 3) {
throw new IllegalArgumentException(
"Illegal plot specification "
+ paramlist
+ "\nPlots must be specified as Q|File|Title"
+ "\nMultiple plots can be separated by ';'.");
}
PlotParameters pp = new PlotParameters();
pp.x = Enum.valueOf(Quantity.class, params[0]);
pp.output = new File(params[1]);
if (params.length == 3) {
pp.title = params[2];
}
lpp.add(pp);
}
kmls = lpp.toArray(new PlotParameters[lpp.size()]);
}
/* text option */
opt = OPT_TEXT;
c = opt.charAt(0);
if (cl.hasOption(c)) {
app.enableOutput(Output.Text, cl.getOptionValue(c));
}
/* no-stats option */
opt = OPT_NO_STATS;
c = opt.charAt(0);
if (cl.hasOption(c)) {
app.setPrintStats(false);
}
/* US units option */
opt = OPT_US_UNITS;
c = opt.charAt(0);
if (cl.hasOption(c)) {
app.setUseUSUnits(true);
}
/* speedcutoff option */
opt = OPT_MIN_SPEED;
c = opt.charAt(0);
if (cl.hasOption(c)) {
app.setSpeedCutoff(Float.parseFloat(cl.getOptionValue(c)));
}
return cl;
}
}
private static class BearingFunction implements UnivariateRealFunction {
private final UnivariateRealFunction dLon;
private final UnivariateRealFunction lat;
private final UnivariateRealFunction dLat;
public BearingFunction(UnivariateRealFunction lonSpeed,
UnivariateRealFunction latSpeed, UnivariateRealFunction lat) {
dLon = lonSpeed;
dLat = latSpeed;
this.lat = lat;
}
/**
* Calculates bearing with respect to true north.
*
* <pre>
* Latitude'
* bearing = 360° - arctan -------------------------
* Longitude' * cos(Latitude)
* where f' = df/dt
* </pre>
*
* @param t
* the dateTime for wich to calculate the bearing expressed
* as a java.util.Date.getTime() value
* @return the bearing at time t
*/
@Override
public double value(double t) throws FunctionEvaluationException {
double dx = dLon.value(t) * Math.cos(lat.value(t) * Math.PI / 180.);
double dy = dLat.value(t);
double result = Math.atan2(dy, dx);
result = ((2.5 - result / Math.PI) % 2.0) * 180.0;
return result;
}
}
private static double[] linspace(final double xmin, final double xmax,
final int n) {
double[] x = new double[n];
final double xstep = (xmax - xmin) / (n - 1);
for (int i = 0; i < n; i++) {
x[i] = xmin + i * xstep;
}
return x;
}
}