BeginPackage["GPXTrack`", {"Vincenty`", "WebServices`"}];

GetTracks::usage = "GetTracks[gpxfile] gets the track XMLElements stored in an XML GPX gpxfile."

GetTrackData::usage = "GetTrack[track] returns a list with elements {lat, lon, elev, time} for each point in track stored as a XMLElement."

AnalyzeTrack::usage = "AnalyzeTrack[track] gives some summary data for a track stored as a XMLElement."

Options[AnalyzeTrack] = {Units->"English", RestThreshold->Automatic};

Options[GetTrackData] = {Time->"Seconds"};

Begin["`Private`"];

InstallService["http://terraservice.net/TerraService.asmx?WSDL"];

$theme = "Topo";
$scale = "Scale16m";

GetTracks[file_] :=
Module[{xml, track},
    xml = Import[file, "XML"];
    track = Cases[xml, XMLElement["trk", ___], Infinity]
];

GetTrackData[track:XMLElement["trk", ___], opts___] := 
Module[{lat, lon, elev, times, get, tform},
    get[pat_] := Cases[track, pat, Infinity];
    lat = get[XMLElement["trkpt", {___, "lat" -> l_, ___ }, __] :> ToExpression[l]];
    lon =  get[XMLElement["trkpt", {___, "lon" -> l_, ___}, __] :> ToExpression[l]];
    elev = get[XMLElement["ele", _, {meters_}] :> ToExpression[meters]];
    times = get[XMLElement["time", {}, {t_}]:> t];
    tform = Time /. Flatten[{opts, Options[GetTrackData]}];
    If[tform =!= "Full",
        times = Map[FromDate[DateFromString[#]]&, times];
        times -= First[times];
        (* To correct for Magellan or GPSExpert bug *)
        times = Map[If[# < 0, Module[{t = #, day = 3600*24}, While[t < 0, t += day];t],#]&, times];
        Switch[tform,
            "Hours", times /= 3600,
            "Minutes", times /= 60];
    ];
    Transpose[{lat, lon, elev, times}]
];

AnalyzeTrack[track:XMLElement["trk", ___], opts___] := 
Module[{data = GetTrackData[track, Time->"Seconds"]},
    {u, reth} = {Units, RestThreshold} /. Flatten[{opts, Options[AnalyzeTrack]}];
    If[u =!= "Metric", u = "English"];
    If[reth === Automatic, 
        reth = If[u === "English", 1., 1.60935]];
    coords = data[[All, {1,2}]];
    showTrack[data[[All,{2,1}]]];
    slengths = First[VincentyDistance[Transpose[Drop[coords, -1]], Transpose[Drop[coords, 1]]]];
    slengths = ConvertDistance[slengths, u];
    dist = FoldList[Plus, 0., slengths];
    Print["Total distance ", dunits[Last[dist], u]];
    h = ConvertElevation[data[[All, 3]], u];
    dh = Drop[h, 1] - Drop[h, -1];
    Print["Net elevation change ", hunits[Total[dh], u]];
    Print["Total climb ", hunits[Total[Select[dh, Positive]], u]];
    Print["Total descent ", hunits[Total[Select[dh, Negative]], u]];
    ListPlot[Transpose[{dist, h}], Frame->True, FrameLabel->{dunits["distance", u], hunits["elev", u], "Elevation Profile", None}];
    times = data[[All,4]];
    Print["Total elapsed time ", showtime[Last[times]]];
    Print["Average speed ", showspeed[3600 Last[dist]/Last[times], u]];
    intervals = Drop[times, 1] - Drop[times, -1];
    speed = 3600 slengths/intervals;
    moving = Position[speed, x_ /; x > reth];
    mint = Extract[intervals, moving];
    mtime = Total[mint];
    Print["Estimated moving time ", showtime[mtime]];
    Print["Estimated average moving speed ", showspeed[3600 Last[dist]/mtime, u]];
    
];
    

DateFromString[sdate_String] := 
Module[{pos},
    pos = StringPosition[sdate, {"-", "T", ":", "Z"}];
    seg[{ps1_, ps2_}, {pe1_, pe2_}] := 
        ToExpression[StringTake[sdate, {ps2 + 1, pe1 - 1}]];
    MapThread[seg, {Join[{{0,0}}, Drop[pos, -1]], pos}]
];


MetersToFeet[m_] := m (100/(12 2.54));
MetersToMiles[m_] := m (100/(12 2.54 5280));
MilesToMeters[m_] := m/(100/(12 2.54 5280));

FromMetersPerSecond[s_, "Metric"] := 3.6 s;
FromMetersPerSecond[s_, "English"] := 3600 MetersToMiles[s];

ConvertElevation[h_, "Metric"] := h;
ConvertElevation[h_, "English"] := MetersToFeet[h];

hunits[s_String, "Metric"] := s<>" (m)";
hunits[s_?NumberQ, "Metric"] := ToString[s]<>" m";
hunits[s_String, "English"] := s<>" (feet)";
hunits[s_?NumberQ, "English"] := ToString[s]<>" feet";
hunits[s_, u_] := hunits[ToString[s], u];

ConvertDistance[d_, "Metric"] := d/1000;
ConvertDistance[d_, "English"] := MetersToMiles[d];

dunits[s_String, "Metric"] := s<>" (km)";
dunits[s_?NumberQ, "Metric"] := ToString[s]<>" km";
dunits[s_String, "English"] := s<>" (mi)";
dunits[s_?NumberQ, "English"] := ToString[s]<>" miles";
dunits[s_, u_] := dunits[ToString[s], u];


showspeed[speed_, "Metric"] := ToString[speed]<>" km/h";
showspeed[speed_, "English"] := ToString[speed]<>" mi/h";

lz[ms_ /; ms < 10] := "0"<>ToString[ms];
lz[ms_] := ToString[ms];

showtime[t_] := Module[{tint, h, m, s},
    tint = Round[t];
    h = Floor[tint/3600];
    m = Floor[(tint - 3600 h)/60];
    s = (tint - 3600 h - 60 m);
    ToString[h]<>":"<>lz[m]<>":"<>lz[s]
];

showTrack[lonlat_] :=
  Module[{llon, llat, nwlon, selon, nwlat, selat, rect, nw, se, nwx, nwy, 
      sex, sey, scene}, 
    llon = lonlat[[All, 1]];
    llat = lonlat[[All, 2]];
    nwlon = Min[llon];
    selon = Max[llon];
    nwlat = Max[llat];
    selat = Min[llat];
    rect = 
      GetAreaFromRect[{"Lon" -> nwlon, "Lat" -> nwlat}, {"Lon" -> selon, 
          "Lat" -> selat}, $theme, $scale];
    nw = "NorthWest" /. rect;
    se = "SouthEast" /. rect;
    {nwx, nwy} = {"X", "Y"} /. ("Id" /. ("TileMeta" /. nw));
    {sex, sey} = {"X", "Y"} /. ("Id" /. ("TileMeta" /. se));
    scene = "Scene" /. ("Id" /. ("TileMeta" /. nw));
    nt = Max[Abs[nwy - sey], Abs[nwx - sex]];
    images = Outer[
        ImportString[
            getImage[{"Theme" -> $theme, "Scale" -> $scale, "Scene" -> scene, 
                "Y" -> #1, "X" -> #2}], "JPEG"] &, Reverse[Range[sey, nwy]], 
        Range[nwx, sex]];
    meta = 
      Outer[GetTileMetaFromTileId[{"Theme" -> $theme, "Scale" -> $scale, 
              "Scene" -> scene, "Y" -> #1, "X" -> #2}] &, 
        Reverse[Range[sey, nwy]], Range[nwx, sex]];
    points = 
      Map[scalePoints[
              Transpose[{"Lon", "Lat"} /. ({"SouthWest", "NorthEast"} /. #)]][
            lonlat] &, meta, {2}];
    both = 
      MapThread[
        Block[{$DisplayFunction = Identity}, Show[#1, #2]] &, {images, 
          points}, 2];
    Show[GraphicsArray[both], GraphicsSpacing -> {0, 0}]
    ];

getImage[id_List] :=
    getImage[id] =
      Block[{
          native = InvokeServiceOperation::native, 
          rspns = InvokeServiceOperation::rspns,
          result},
        Off[InvokeServiceOperation::native];
        Off[InvokeServiceOperation::rspns];
        result = GetTile[id];
        If[result === $Failed, result = $Failed];
        If[Head[native] =!= $Off, On[InvokeServiceOperation::native]];
        If[Head[rspns] =!= $Off, On[InvokeServiceOperation::rspns]];
        result
        ];

makeCond = 
    Function[Condition[{plon_, plat_}, 
        Evaluate[
          And[ #[[1, 1]] <= plon <= #[[1, 2]], #[[2, 1]] <= 
              plat <= #[[2, 2]]]]]];

scalePoints[{{lon1_, lon2_}, {lat1_, lat2_}}][data_] := 
  Block[{selected, dlon, dlat, plon, plat},
    selected = Cases[data, makeCond[{{lon1, lon2}, {lat1, lat2}}]];
    If[selected === {}, 
      Graphics[selected],
      {dlon, dlat} = Transpose[selected];
      dlon = Rescale[dlon, {lon1, lon2}, {0, 199}];
      dlat = Rescale[dlat, {lat1, lat2}, {1, 200}];
      Graphics[{RGBColor[1, 0, 0], PointSize[0.01 nt], 
          Map[Point, Transpose[{dlon, dlat}]]}]]];

    
End[];

EndPackage[];
