I am using this code to python and at the end I have calculated the velocity and now I want to take the control of velocity that it should not be bigger then threshold and should do new matching please help me in this regard
Code:
# -*- coding: cp1252 -*-
# Import native arcgisscripting module
import sys, math,arcpy
#from pylab import *
def match_coo(XA,YA,XE,YE,x_gps,y_gps):
global secmin
global Y_end # Retransformed, matched coordinates
global X_end
global x_min # minimal distance to the edge
beta=-math.atan((XE-XA)/(YE-YA))
o=math.sin(beta)
a=math.cos(beta)
x0=-o*YA-a*XA
y0=-a*YA+o*XA
xi=x0+o*y_gps+a*x_gps
yi=y0+a*y_gps-o*x_gps
dis=((XE-XA)**2+(YE-YA)**2)**0.5
if abs(xi)<x_min:
if yi<dis and yi>0:
secmin=x_min
x_min=abs(xi)
xi=0
Y_end=((a)/(a**2+o**2))*(yi-y0)+((o)/(a**2+o**2))*(xi-x0)
X_end=((-o)/(a**2+o**2))*(yi-y0)+((a)/(a**2+o**2))*(xi-x0)
if abs(xi)<secmin and abs(xi)>x_min:
secmin=xi
def writing_fc(x,y,cursor):
feat = cursor.newRow()
pnt.X=x
pnt.Y=y
# Set the geometry of the new feature to the array of points
feat.shape = pnt
cursor.insertRow(feat)
######################################################################
######################################################################
#MAIN program
#arcpy.OverWriteOutput = True
arcpy.env.overwriteOutput = True
arcpy.AddMessage('----------------------------------')
arcpy.AddMessage('IIGS, University of Stuttgart 2011')
arcpy.AddMessage('----------------------------------')
arcpy.AddMessage('starting Map Matching Script...')
try:
paramGDB=str(sys.argv[1])
paramInStreets=str(sys.argv[2])
paramIngps=str(sys.argv[3])
paramBuffer=50
temp=str(paramGDB+ "")
except Exception, ErrorDesc:
arcpy.AddError("Input parameters are incorrect")
arcpy.env.Workspace=paramGDB
a= 9999
o= 9999
y_min=9999
x_min=9999
secmin=9999
X_end=0
Y_end=0
X_end_0=0
Y_end_0=0
i=0
threshold=3
arcpy.MakeFeatureLayer_management(paramInStreets, "streets_lyr")
arcpy.MakeFeatureLayer_management(paramIngps, "gps_lyr")
gps_id=[]
gps_x=[]
gps_y=[]
time=[]
vx=[]
vy=[]
vx1=[]
vy1=[]
count_gps=0
# read gps points
# --------------------------------------
desc_g = arcpy.Describe(paramIngps)
shapefieldname_g = desc_g.ShapeFieldName
fielder=desc_g.fields
rows_gps = arcpy.SearchCursor(paramIngps)
for row in rows_gps:
# Create the geometry object 'feat'
#
feat2 = row.getValue(shapefieldname_g)
pnt = feat2.getPart()
t = row.time
# Print x,y coordinates of current gps point
#
arcpy.AddMessage("gps: "+ str(pnt.X)+" "+str(pnt.Y))
gps_id.append(pnt.ID)
gps_x.append(pnt.X)
gps_y.append(pnt.Y)
time.append(t)
#time.diff(t)
del rows_gps
count_gps=len(gps_id)
print count_gps
arcpy.AddMessage('1. --> '+str(count_gps)+' gps Points from "'+paramIngps+'" are imported')
arcpy.AddMessage("--------------------------------------------------------")
#select only the features of streets Layer which are closer than 50 m to the gps points
arcpy.SelectLayerByLocation_management("streets_lyr", "WITHIN_A_DISTANCE", "gps_lyr", paramBuffer)
#copy the selected features to a new feature class
arcpy.CopyFeatures_management("streets_lyr",temp)
#arcpy.MakeFeatureLayer_management(temp,"tmp_lyr")
arcpy.AddMessage('2. New Layer with selected streets elements near the gps points is created. Distance='+str(paramBuffer)+' m\n')
#template = arcpy.GetParameterAsText(2)
# new Feature Class
fcname="\matchedPoi.shp"
arcpy.AddMessage(paramGDB+" "+fcname+" creating...")
sr = arcpy.Describe(paramIngps).spatialReference
arcpy.CreateFeatureclass_management(paramGDB,fcname, "POINT", "", "DISABLED", "DISABLED",sr)
fcname=paramGDB+"/"+fcname
arcpy.MakeFeatureLayer_management(fcname, "results")
arcpy.AddMessage("Done!")
cur = arcpy.InsertCursor(fcname)
pnt = arcpy.CreateObject("Point")
#read Streets
#----------------------------------------
# Identify the geometry field
desc = arcpy.Describe("streets_lyr")
shapefieldname = desc.ShapeFieldName
arcpy.AddMessage('3. reading streets elements...\n')
for i in range(0,count_gps):
rows = arcpy.SearchCursor("streets_lyr")
for row in rows:
# Create the geometry object
feat = row.getValue(shapefieldname)
partnum = 0
# Count the number of points in the current multipart feature
partcount = feat.partCount
# Enter while loop for each part in the feature (if a singlepart feature
# this will occur only once)
while partnum < partcount:
# Print the part number
#print "Part " + str(partnum) + ":"
part = feat.getPart(partnum)
pntcount = 0
# Enter while loop for each vertex
for pnts in part:
if pntcount!=0:
#arcpy.AddMessage("Function --> ")
#here the function trafo_par is executed !!!!!!
# --------------------------------------------
match_coo(x0,y0,pnts.X,pnts.Y,gps_x[i],gps_y[i])
# ---------------------------------------------
#arcpy.AddMessage("Function end! ")
# If pnt is null, either the part is finished or there is an
# interior ring
x0=pnts.X
y0=pnts.Y
pntcount += 1
if not pnts:
print ''
if pnts:
print "Interior Ring:"
partnum += 1
arcpy.AddMessage("--> matched point: "+ str(X_end) +" ; "+ str(Y_end))
arcpy.AddMessage("--> time : "+ str(time[i]))
vx=(X_end-X_end_0)/5
vy=(Y_end-Y_end_0)/5
X_end_0=X_end
Y_end_0=Y_end
arcpy.AddMessage("velocity for each matched point "+ str(vx) +" ; "+ str(vy))
#plot(vx,vy)
#show(plot)
writing_fc(X_end,Y_end,cur)
#writing_vel(vx,vy)
x_min=99999
arcpy.AddMessage("--------------------------------------------------------")
arcpy.AddMessage('All points have been matched ')
arcpy.Delete_management(temp)
for vx in range(vx):
if vx<threshold:
vx1=vx
return vx
for vy in range(vy):
if vy<threshold:
vy1=vy
return vy
arcpy.AddMessage("--> newmatched point: "+ str(X_end) +" ; "+ str(Y_end))