Initial commit

gpxanalyse.py

+import xml.dom.minidom
+from matplotlib.pyplot import plot, show
+from pylab import polyfit, polyval
+from tkFileDialog import askopenfilename
+from math import pi, sin, sqrt, cos, atan, fabs
+
+def dms(time):
+    s=int(time%60)
+    m=int(time/60%60)
+    h=int(time/3600)
+    return (h,m,s)
+
+def dist(lat1, lon1, lat2, lon2):
+    if (lat1,lon1)==(lat2,lon2): return 0
+    f = 1.0/298.2572221
+    a = 6378137.0
+    F = (lat1+lat2)/360.0*pi
+    G = (lat1-lat2)/360.0*pi
+    l = (lon1-lon2)/360.0*pi
+    S = sin(G)**2 * cos(l)**2 + cos(F)**2 * sin(l)**2
+    C = cos(G)**2 * cos(l)**2 + sin(F)**2 * sin(l)**2
+    w = atan(sqrt(S/C))
+    R = sqrt(S*C)/w
+    D = 2.0*w*a
+    H1 = (3*R-1)/(2*C)
+    H2 = (3*R+1)/(2*S)
+    return D*(1.0 + f*H1*sin(F)**2 * cos(G)**2 - f*H2*cos(F)**2 * sin(G)**2)
+
+def interpol(x, y, lx):
+    i=0
+    while lx > x[i]:
+        i+=1
+
+    lu = x[i-1]
+    lo = x[i]
+    hu = y[i-1]
+    ho = y[i]
+
+    if lu==lo:
+        hx = hu
+    else:
+        hx = hu+ (ho-hu)/(lo-lu)*(lx-lu)
+    return hx
+
+def handleGpx(gpx):
+    class Data:
+        lat = []
+        lon = []
+        tim = []
+        ele = []
+        dis = []
+
+    s=0.0
+    i=0
+    pts=Data();
+    tps = gpx.getElementsByTagName("trkpt")
+    for tp in tps:
+        tpv = handleTrackpoint(tp)
+        if tpv:
+            (tim,ele,lat,lon) = tpv
+            pts.tim.append(tim)
+            pts.ele.append(ele)
+            pts.lat.append(lat)
+            pts.lon.append(lon)
+            if i>0:
+                pts.dis.append(pts.dis[i-1]+dist(pts.lat[i-1],pts.lon[i-1],pts.lat[i],pts.lon[i]))
+            else:
+                pts.dis.append(0)
+            i=i+1
+    return pts
+
+def handleTrackpoint(tp):
+    lat = float(tp.getAttribute("lat"))
+    lon = float(tp.getAttribute("lon"))
+    tim=0
+    for t in tp.getElementsByTagName("time")[0].childNodes[0].data.split('T')[1].split(':'):
+        tim = tim*60 + int(t[:2])
+    ele = (float(tp.getElementsByTagName("ele")[0].childNodes[0].data))
+    return tim, ele, lat, lon
+
+def dpeuker(x,y,eps):
+    xs = []
+    ys = []
+    xs.append(x[0])
+    ys.append(y[0])
+    dp(x,y,eps,xs,ys,0,len(x)-1)
+    xs.append(x[-1])
+    ys.append(y[-1])
+    return xs,ys
+
+def dp(x,y,eps,xs,ys,iu,io):
+    if iu+1 < io:
+        maxi=0
+        maxd=0.0
+        for i in range(iu+1,io):
+            dx = x[io] - x[iu]
+            dy = y[io] - y[iu]
+            d = fabs(y[i] - y[iu] - dy/dx * (x[i]-x[iu]))
+            if d>maxd:
+                maxi=i
+                maxd=d
+        if maxd>eps:
+            dp(x,y,eps,xs,ys,iu,maxi)
+            xs.append(x[maxi])
+            ys.append(y[maxi])
+            dp(x,y,eps,xs,ys,maxi,io)
+
+
+filename = askopenfilename()
+if len(filename)>4: 
+    print filename
+    dom = xml.dom.minidom.parse(filename)
+    data = handleGpx(dom)
+
+    plot(data.lon,data.lat)
+    show()
+
+    x,y = dpeuker(data.dis,data.ele,3.0)
+
+    plot(data.dis,data.ele,'b.')
+    plot(x,y,'r-')
+    show()
+

simplify.py

+import xml.dom.minidom
+from matplotlib.pyplot import plot, show
+from pylab import polyfit, polyval
+from tkFileDialog import askopenfilename
+from math import pi, sin, sqrt, cos, atan, fabs
+
+class Point:
+    def __init__(self, latitude, longitude, time, elevation, speed=0, distance=0, segmentlength=0):
+        self.latitude = latitude
+        self.longitude = longitude
+        self.time = time
+        self.elevation = elevation
+        self.distance = distance
+        self.speed = speed
+        self.segmentlength = segmentlength
+
+    @staticmethod
+    def distanceBetween(p1, p2):
+        if (p1.latitude,p1.longitude)==(p2.latitude,p2.longitude): return 0
+        f = 1.0/298.2572221
+        a = 6378137.0
+        F = (p1.latitude+p2.latitude)/360.0*pi
+        G = (p1.latitude-p2.latitude)/360.0*pi
+        l = (p1.longitude-p2.longitude)/360.0*pi
+        S = sin(G)**2 * cos(l)**2 + cos(F)**2 * sin(l)**2
+        C = cos(G)**2 * cos(l)**2 + sin(F)**2 * sin(l)**2
+        w = atan(sqrt(S/C))
+        R = sqrt(S*C)/w
+        D = 2.0*w*a
+        H1 = (3*R-1)/(2*C)
+        H2 = (3*R+1)/(2*S)
+        return D*(1.0 + f*H1*sin(F)**2 * cos(G)**2 - f*H2*cos(F)**2 * sin(G)**2)
+
+    def distanceTo(self, point):
+        return Point.distanceBetween(self, point)
+
+# def dms(time):
+#     s=int(time%60)
+#     m=int(time/60%60)
+#     h=int(time/3600)
+#     return (h,m,s)
+
+def interpol(x, y, lx):
+    i=0
+    while lx > x[i]:
+        i+=1
+
+    lu = x[i-1]
+    lo = x[i]
+    hu = y[i-1]
+    ho = y[i]
+
+    if lu==lo:
+        hx = hu
+    else:
+        hx = hu+ (ho-hu)/(lo-lu)*(lx-lu)
+    return hx
+
+def handleGpx(gpx):
+    points = [];
+    tps = gpx.getElementsByTagName("trkpt")
+    for i in xrange(len(tps)):
+        tp = tps[i]
+        (tim,ele,lat,lon) = handleTrackpoint(tp)
+        point = Point(latitude=lat, longitude=lon, elevation=ele, time=tim)
+        if i>0:
+            point.segmentlength = point.distanceTo(points[i-1])
+            point.distance = points[i-1].distance + point.segmentlength
+            point.speed = point.segmentlength / (point.time - points[i-1].time)
+        points.append(point)
+    return points
+
+def handleTrackpoint(tp):
+    lat = float(tp.getAttribute("lat"))
+    lon = float(tp.getAttribute("lon"))
+    tim=0
+    for t in tp.getElementsByTagName("time")[0].childNodes[0].data.split('T')[1].split(':'):
+        tim = tim*60 + int(t[:2])
+    ele = (float(tp.getElementsByTagName("ele")[0].childNodes[0].data))
+    return tim, ele, lat, lon
+
+def dpeuker(x,y,eps):
+    xs = []
+    ys = []
+    xs.append(x[0])
+    ys.append(y[0])
+    dp(x,y,eps,xs,ys,0,len(x)-1)
+    xs.append(x[-1])
+    ys.append(y[-1])
+    return xs,ys
+
+def dp(x,y,eps,xs,ys,iu,io):
+    if iu+1 < io:
+        maxi=0
+        maxd=0.0
+        for i in range(iu+1,io):
+            dx = x[io] - x[iu]
+            dy = y[io] - y[iu]
+            d = fabs(y[i] - y[iu] - dy/dx * (x[i]-x[iu]))
+            if d>maxd:
+                maxi=i
+                maxd=d
+        if maxd>eps:
+            dp(x,y,eps,xs,ys,iu,maxi)
+            xs.append(x[maxi])
+            ys.append(y[maxi])
+            dp(x,y,eps,xs,ys,maxi,io)
+
+
+filename = "ol_egestorf.gpx"
+dom = xml.dom.minidom.parse(filename)
+data = handleGpx(dom)
+
+for p in data:
+    print p.segmentlength, p.speed