|
Salt Lake City uses GPS to help gather travel time data to help with traffic issues. After converting from .gpx to a feature class, I have the X, Y, and timestamp on each point from the track log. The code below will add fields that are populated with the number of seconds between points, as well as the speed between a point and it's preceeding point. This code is kind of fun. Calling the addGPSdata function on the feature class first checks for the needed additional data fields and adds them if needed. Next a dictionary of FID, X, Y, Time is built up in memory for accessing later. A cursor is created to update all of the rows in the feature class and for each row, the time and distance is compared and written into the appropriate field. NOTE: all of the time field data's parsing is hardcoded and would need to be altered for a specific time format, and the cursors are hardcoded as well I duct taped a GPS to my dog and he can run 42 MPH!
######################################################### #Name: deriveSpeedFromGPS.py #Author: Kevin Bell #Date: 20080108 #Purpose: add speed/seconds to gps data import time import arcgisscripting #-------------------------------------- def _getTimeCode(inFID,d): # TIP: google for python string slicing! '''Parse a time & return seconds from epoch''' t = d[inFID][2] y = int(t[:4]) m = int(t[5:7]) day = int(t[8:10]) h = int(t[11:13]) minute = int(t[14:16]) sec = int(t[17:19]) timestamp = y,m,day,h,minute,sec,-1, -1, -1 timeCode = time.mktime(timestamp) return timeCode def _getDistanceFromPreceeding(inFID,d): # TIP: google for python dictionaries '''return the distance from the point that preceeds the current point''' x1 = d[inFID][0] y1 = d[inFID][1] inFID2 = inFID - 1 x2 = d[inFID2][0] y2 = d[inFID2][1] xDelta = abs(x2 - x1) yDelta = abs(y2 - y1) return round(xDelta + yDelta, 2) def _buildXYTdict(fc): # TIP: see the GP cursor help topics '''Build a dict of XYT's needed later to compare point against point''' d = {} cur = gp.SearchCursor(fc) row = cur.Next() while row: d[row.FID] = [row.Y_PROJ, row.X_PROJ, row.TIME] row = cur.Next() del cur return d def _checkFields(fc, field): '''check for fields needed to hold new data, add them if necessary ''' fldList = gp.ListFields(fc, field) if fldList.Next() == None: gp.AddField_management(fc, field, "FLOAT", 9, 2) print "added field: " + field def addGPSdata(fc): '''add SpeedMPH and seconds to FC... this is relative to its preceeding point''' _checkFields(fc, "seconds") _checkFields(fc, "SpeedMPH") d = _buildXYTdict(fc) updcur = gp.UpdateCursor(fc) updrow = updcur.Next() while updrow: theFID = updrow.FID if theFID != 0: # if it's the 1st point, there's nothing to compare to... y = theFID - 1 curTC = _getTimeCode(theFID, d) prevTC = _getTimeCode(y,d) secs = curTC - prevTC; print str(secs) + " seconds" meters = _getDistanceFromPreceeding(theFID, d); print str(meters) + " meters" speed = (meters/secs) * 2.23693629; print str(round(speed,2)) + " mph" + "\n"
updrow.seconds = secs updrow.SpeedMPH = round(speed,2) updcur.UpdateRow(updrow) updrow = updcur.Next() else: updrow = updcur.Next() del updcur #---------------------------------------------------------- gp = arcgisscripting.create() gp.workspace = r"C:\dog" fc = "dogNAD83_UTM12N.shp"
addGPSdata(fc); print "...added GPS data" del gp print "...all done." |