Expand|Select|Wrap|Line Numbers
- # -*- 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))