By using this site, you agree to our updated Privacy Policy and our Terms of Use. Manage your Cookies Settings.
464,728 Members | 1,145 Online
Bytes IT Community
+ Ask a Question
Need help? Post your question and get tips & solutions from a community of 464,728 IT Pros & Developers. It's quick & easy.

Map Matching script

P: 1
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
Expand|Select|Wrap|Line Numbers
  1. # -*- coding: cp1252 -*-
  2. # Import native arcgisscripting module
  3. import  sys, math,arcpy
  4. #from pylab import *
  5.  
  6.  
  7. def match_coo(XA,YA,XE,YE,x_gps,y_gps):
  8.  
  9.     global secmin
  10.     global Y_end      # Retransformed, matched coordinates
  11.     global X_end
  12.     global x_min      # minimal distance to the edge
  13.     beta=-math.atan((XE-XA)/(YE-YA))
  14.     o=math.sin(beta)
  15.     a=math.cos(beta)
  16.     x0=-o*YA-a*XA
  17.     y0=-a*YA+o*XA
  18.     xi=x0+o*y_gps+a*x_gps
  19.     yi=y0+a*y_gps-o*x_gps
  20.     dis=((XE-XA)**2+(YE-YA)**2)**0.5
  21.     if abs(xi)<x_min:
  22.         if yi<dis and yi>0:
  23.             secmin=x_min
  24.             x_min=abs(xi)
  25.             xi=0
  26.             Y_end=((a)/(a**2+o**2))*(yi-y0)+((o)/(a**2+o**2))*(xi-x0)
  27.             X_end=((-o)/(a**2+o**2))*(yi-y0)+((a)/(a**2+o**2))*(xi-x0)
  28.     if abs(xi)<secmin and abs(xi)>x_min:
  29.         secmin=xi
  30.  
  31.  
  32. def writing_fc(x,y,cursor):
  33.     feat = cursor.newRow()
  34.  
  35.     pnt.X=x
  36.     pnt.Y=y
  37.  
  38.     # Set the geometry of the new feature to the array of points
  39.     feat.shape = pnt
  40.     cursor.insertRow(feat)
  41. ######################################################################
  42. ######################################################################
  43. #MAIN program
  44. #arcpy.OverWriteOutput = True
  45. arcpy.env.overwriteOutput = True
  46.  
  47. arcpy.AddMessage('----------------------------------')
  48. arcpy.AddMessage('IIGS, University of Stuttgart 2011')
  49. arcpy.AddMessage('----------------------------------')
  50. arcpy.AddMessage('starting Map Matching Script...')
  51.  
  52. try:
  53.     paramGDB=str(sys.argv[1])
  54.     paramInStreets=str(sys.argv[2])
  55.     paramIngps=str(sys.argv[3])
  56.     paramBuffer=50
  57.     temp=str(paramGDB+ "")
  58.  
  59.  
  60.  
  61. except Exception, ErrorDesc:
  62.     arcpy.AddError("Input parameters are incorrect")
  63.  
  64. arcpy.env.Workspace=paramGDB
  65. a= 9999
  66. o= 9999
  67. y_min=9999
  68. x_min=9999
  69. secmin=9999
  70. X_end=0
  71. Y_end=0
  72. X_end_0=0
  73. Y_end_0=0
  74. i=0
  75. threshold=3
  76.  
  77. arcpy.MakeFeatureLayer_management(paramInStreets, "streets_lyr")
  78. arcpy.MakeFeatureLayer_management(paramIngps, "gps_lyr")
  79.  
  80.  
  81.  
  82. gps_id=[]
  83. gps_x=[]
  84. gps_y=[]
  85. time=[]
  86. vx=[]
  87. vy=[]
  88. vx1=[]
  89. vy1=[]
  90. count_gps=0
  91.  
  92.  
  93.  
  94. # read gps points
  95. # --------------------------------------
  96. desc_g = arcpy.Describe(paramIngps)
  97. shapefieldname_g = desc_g.ShapeFieldName
  98. fielder=desc_g.fields
  99.  
  100. rows_gps = arcpy.SearchCursor(paramIngps)
  101.  
  102.  
  103. for row in rows_gps:
  104.  
  105.      # Create the geometry object 'feat'
  106.      #
  107.      feat2 = row.getValue(shapefieldname_g)
  108.      pnt = feat2.getPart()
  109.      t = row.time
  110.  
  111.  
  112.     # Print x,y coordinates of current gps point
  113.     #
  114.      arcpy.AddMessage("gps: "+ str(pnt.X)+" "+str(pnt.Y))
  115.      gps_id.append(pnt.ID)
  116.      gps_x.append(pnt.X)
  117.      gps_y.append(pnt.Y)
  118.      time.append(t)
  119.      #time.diff(t)
  120.  
  121.  
  122. del rows_gps    
  123. count_gps=len(gps_id)
  124. print count_gps
  125. arcpy.AddMessage('1. --> '+str(count_gps)+' gps Points from "'+paramIngps+'" are imported')
  126. arcpy.AddMessage("--------------------------------------------------------")
  127.  
  128.  
  129. #select only the features of streets Layer which are closer than 50 m to the gps points
  130. arcpy.SelectLayerByLocation_management("streets_lyr", "WITHIN_A_DISTANCE", "gps_lyr", paramBuffer)
  131. #copy the selected features to a new feature class
  132. arcpy.CopyFeatures_management("streets_lyr",temp)
  133. #arcpy.MakeFeatureLayer_management(temp,"tmp_lyr")
  134. arcpy.AddMessage('2. New Layer with selected streets elements near the gps points is created. Distance='+str(paramBuffer)+' m\n')
  135. #template = arcpy.GetParameterAsText(2)
  136.  
  137.  
  138. # new Feature Class
  139. fcname="\matchedPoi.shp"
  140. arcpy.AddMessage(paramGDB+" "+fcname+" creating...")
  141. sr = arcpy.Describe(paramIngps).spatialReference
  142. arcpy.CreateFeatureclass_management(paramGDB,fcname, "POINT", "", "DISABLED", "DISABLED",sr)
  143. fcname=paramGDB+"/"+fcname
  144. arcpy.MakeFeatureLayer_management(fcname, "results")
  145. arcpy.AddMessage("Done!")
  146. cur = arcpy.InsertCursor(fcname)
  147. pnt = arcpy.CreateObject("Point")
  148.  
  149.  
  150.  
  151. #read Streets
  152. #----------------------------------------
  153.  
  154.  
  155. # Identify the geometry field
  156. desc = arcpy.Describe("streets_lyr")
  157. shapefieldname = desc.ShapeFieldName
  158. arcpy.AddMessage('3. reading streets elements...\n')
  159. for i in range(0,count_gps):
  160.     rows = arcpy.SearchCursor("streets_lyr")
  161.     for row in rows:
  162.        # Create the geometry object
  163.        feat = row.getValue(shapefieldname)
  164.        partnum = 0
  165.        # Count the number of points in the current multipart feature
  166.        partcount = feat.partCount
  167.        # Enter while loop for each part in the feature (if a singlepart feature
  168.        # this will occur only once)
  169.        while partnum < partcount:
  170.            # Print the part number
  171.            #print "Part " + str(partnum) + ":"
  172.            part = feat.getPart(partnum)
  173.            pntcount = 0
  174.            # Enter while loop for each vertex
  175.  
  176.            for pnts in part:
  177.                    if pntcount!=0:
  178.                       #arcpy.AddMessage("Function -->  ")
  179.                       #here the function trafo_par is executed !!!!!!
  180.                       # -------------------------------------------- 
  181.                       match_coo(x0,y0,pnts.X,pnts.Y,gps_x[i],gps_y[i])
  182.                       # ---------------------------------------------
  183.                       #arcpy.AddMessage("Function end!  ")
  184.                    # If pnt is null, either the part is finished or there is an 
  185.                    # interior ring
  186.                    x0=pnts.X
  187.                    y0=pnts.Y
  188.                    pntcount += 1
  189.                    if not pnts: 
  190.                      print ''
  191.                      if pnts:
  192.                         print "Interior Ring:"
  193.  
  194.  
  195.  
  196.            partnum += 1
  197.  
  198.     arcpy.AddMessage("--> matched point:  "+ str(X_end) +" ; "+ str(Y_end))
  199.     arcpy.AddMessage("--> time         :  "+ str(time[i]))
  200.     vx=(X_end-X_end_0)/5
  201.     vy=(Y_end-Y_end_0)/5
  202.     X_end_0=X_end
  203.     Y_end_0=Y_end
  204.     arcpy.AddMessage("velocity for each matched point  "+ str(vx) +" ; "+ str(vy))
  205.  
  206.  
  207.  
  208.  
  209.     #plot(vx,vy)
  210.     #show(plot)
  211.  
  212.     writing_fc(X_end,Y_end,cur)
  213.     #writing_vel(vx,vy)
  214.  
  215.  
  216.  
  217.  
  218.  
  219.     x_min=99999
  220. arcpy.AddMessage("--------------------------------------------------------")
  221. arcpy.AddMessage('All points have been matched  ') 
  222. arcpy.Delete_management(temp)
  223. for vx in range(vx):
  224.         if vx<threshold:
  225.             vx1=vx
  226.             return vx
  227. for vy in range(vy):
  228.     if vy<threshold:
  229.             vy1=vy
  230.             return vy
  231.  
  232. arcpy.AddMessage("--> newmatched point:  "+ str(X_end) +" ; "+ str(Y_end))
Mar 28 '12 #1
Share this question for a faster answer!
Share on Google+

Post your reply

Sign in to post your reply or Sign up for a free account.