6 Extra

Following steps were performed to segment the road into 500m sections and to calculate the speed, turtuosity and elevation difference for each section.

  1. Open QGIS and add vector layer (select the gpx file)

  2. Select route_points layer

  3. From ‘Menu->Plugins’ select ‘python console’. Click ‘Show editor’ tool in the console. Copy the following code and paste to the editor. Save file as python (such as:gpx.py) file.

  4. Run python code from console.

Code:

from PyQt4.Qt import QVariant
from PyQt4.QtGui import QColor
def getDistance(p1, p2):
    m = d.measureLine(p1, p2)
    m = d.convertMeasurement(m, 2, 0, False)
    return m[0]
def getLineLength(l):
    m = d.measureLength(l)
    m = d.convertMeasurement(m, 2, 0, False)
    return m[0]
def getRadius(pArr, ll):
    pt = []
    pt.append(pArr[0])
    l = len(pArr)
    pt.append(pArr[l-1])
    pline = QgsGeometry.fromPolyline(pt)
    return ll/getLineLength(pline)
    
def interpolate_pt(p1, p2,  d, d1):
    f = (d - d1)/d
    x1 = p1.x()
    x2 = p2.x()
    y1 = p1.y()
    y2 = p2.y()
    x = x1 + (x2 - x1)*f
    y = y1 + (y2 - y1)*f
    
    # print x1, x2, d, d1, x, y
    return QgsPoint(x, y)
    
# create split line
track = iface.activeLayer()

vLyr = QgsVectorLayer('linestring?crs=EPSG:4326','track100','memory')
QgsMapLayerRegistry.instance().addMapLayer(track)
vpr = vLyr.dataProvider()
vpr.addAttributes([QgsField("pid", QVariant.Int),
                   QgsField("time", QVariant.Double),
                   QgsField("length",QVariant.Double),
                   QgsField("speed",QVariant.Double),
                   QgsField("elevdif",QVariant.Double),
                   QgsField("radius", QVariant.Double), # radius ratio
                   QgsField("fspeed", QVariant.Double),
                   QgsField("ftime", QVariant.Double)]) 
vLyr.updateFields()

feat = track.getFeatures()
trackpt = []
flen = 0.0
for f in feat:
    trackpt.append(f)
npt = len(trackpt)
# starting point (Initial points can be skipped)
fp = 500
# end point (Last few points can be skipped)
ep = 500
#length of each segment
seg = 500.0
p1 = trackpt[fp].geometry().asPoint()
pts = [p1]
j = 0
vF = []
dt = 0
dAlt = 0.0
# coefficients calculated from R analysis
coef = {'elev':-1.049, 'radius':-19.466,'intercept':95.235}
for i in range(fp + 1, npt - ep ):
    d = QgsDistanceArea()
    p2 = trackpt[i   ].geometry().asPoint()
    dist = getDistance(p1, p2)
    flen = flen + dist
    pts.append(p2)
    # time 
    tm =( trackpt[i].attribute("time")).toTime_t() - ( trackpt[i -1].attribute("time")).toTime_t() 
    dt += tm
    dAlt += abs(trackpt[i].attribute("ele") - trackpt[i -1].attribute("ele") )
    if flen > seg:
        dd = flen - seg
        p1 = interpolate_pt( p1, p2, dist, dd)
        pts.pop()
        pts.append(p1)
        dt = dt - tm + (dd / dist * tm)
        j += 1
        pline = QgsGeometry.fromPolyline(pts)
        flen = getLineLength(pline)
        qf = QgsFeature()
        qf.setGeometry(pline)
        speed = flen /dt * 3.6
        radius = getRadius(pts, flen)
        rad2 = radius
        if rad2 > 2.0:
            rad2 = 2.0
        fspeed = coef['elev']*dAlt + coef['radius']*rad2+coef['intercept']
        ftime = 1800.0 / fspeed
        qf.setAttributes([j, dt, flen, speed, dAlt,radius,fspeed, ftime])
        vF.append(qf)
        # print j, i, flen
        flen = 0
        pts  = [p1]
        i = i -1
        dt = 0
        dAlt = 0.0
    else:
        p1 = p2
vpr.addFeatures(vF)
vLyr.commitChanges()
filename = "path/to/tracksp.shp"
QgsVectorFileWriter.writeAsVectorFormat(vLyr, filename,
                "UTF-8", None, "ESRI Shapefile")

vLyrB =  QgsVectorLayer(filename, "Speed(actual)", "ogr")

vLyr = QgsVectorLayer(filename, "GPX track (speed)", "ogr")
# Render
speedrange = (
("Very Low (0 ~ 10km/hr)", 0.0, 10.0, "cyan"),
("Low(10 ~ 20km/hr)", 10.0001, 20.0, "yellow"),
("Moderate (20 ~ 40km/hr)", 20.0001, 40.0, "green"),
("High (40 ~ 60km/hr)", 40.0001, 60.0, "orange"),
("Very High ( > 60km/hr)", 60.0001, 200.0, "red"))

ranges = []
ranges2 = []
for label, lower, upper, color in speedrange:
    sym = QgsSymbolV2.defaultSymbol(vLyr.geometryType())
    sym.setColor(QColor(color))
    sym.setWidth(0.6)
    rng = QgsRendererRangeV2(lower, upper, sym, label)
    ranges.append(rng)
    sym.setWidth(1.5)
    rng = QgsRendererRangeV2(lower, upper, sym, label)
    ranges2.append(rng)

field = "speed"
renderer = QgsGraduatedSymbolRendererV2(field, ranges)
vLyr.setRendererV2(renderer)
renderer = QgsGraduatedSymbolRendererV2('fspeed', ranges2)
vLyrB.setRendererV2(renderer)
# label

lbl = QgsPalLayerSettings()
lbl.readFromLayer(vLyrB)
lbl.enabled = True
lbl.isExpression = True
lbl.fieldName = "round(\"speed\", 0) || round(\"fspeed\",0)"
lbl.bufferDraw = True
lbl.bufferColor = QColor("blue")
lbl.bufferSize = 0.7
lbl.textColor = QColor("white")
lbl.placement = QgsPalLayerSettings.AboveLine
lbl.setDataDefinedProperty(QgsPalLayerSettings.Size, True, True, '8', ' ')
lbl.writeToLayer(vLyrB)
QgsMapLayerRegistry.instance().addMapLayer(vLyrB)
QgsMapLayerRegistry.instance().addMapLayer(vLyr)
print "finished"