473,320 Members | 1,914 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,320 software developers and data experts.

Map Matching script

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
0 1630

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

Similar topics

1
by: Dariusz | last post by:
Despite looking at a number of tutorials in books and online examples on pattern matching - I can't quite get my head around it to work... so need some help. What I have now is a variable that...
1
by: Adam | last post by:
Hello all, I have written a script that utilizes mysql_unbuffered_query to quickly perform a series of if statements on the returned row as my result sets are roughly 50,000,000+ in size. ...
2
by: Peter Fein | last post by:
I'm trying to use pyparsing write a screenscraper. I've got some arbitrary HTML text I define as opener & closer. In between is the HTML data I want to extract. However, the data may contain the...
17
by: Andrew McLean | last post by:
I have a problem that is suspect isn't unusual and I'm looking to see if there is any code available to help. I've Googled without success. Basically, I have two databases containing lists of...
0
by: Eric A. Hall | last post by:
I'm trying to scratch up a whois-like client in perl, and am having a strange problem with matching input characters. Specifically, matching against things like "@" works fine over the network, but...
7
by: Niko | last post by:
Hi , I am writing script using multiple pattern matching. I need to replace the $_ =~ few times. Example: sub xx { open (TEMPF,"> $T_FILE");
11
by: rajarao | last post by:
hi I want to remove the content embedded in <script> and </script> tags submitted via text box. My java script should remove the content embedded between <script> and </script> tag. my current...
2
by: Kae Verens | last post by:
Anyone know of a script which returns a list of elements matching a specified CSS selector? Kae
4
by: 28tommy | last post by:
Hi, I'm trying to find scripts in html source of a page retrieved from the web. I'm trying to use the following rule: match = re.compile('<script + src=+>') I'm testing it on a page that...
2
by: jonathan184 | last post by:
Hi I am having a problme where the results of the sql count is not matching the results of the perl script sql count. The script was working fine up till Wed last week and after that the results...
0
by: DolphinDB | last post by:
Tired of spending countless mintues downsampling your data? Look no further! In this article, you’ll learn how to efficiently downsample 6.48 billion high-frequency records to 61 million...
0
by: ryjfgjl | last post by:
ExcelToDatabase: batch import excel into database automatically...
0
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
1
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
0
by: Vimpel783 | last post by:
Hello! Guys, I found this code on the Internet, but I need to modify it a little. It works well, the problem is this: Data is sent from only one cell, in this case B5, but it is necessary that data...
0
by: ArrayDB | last post by:
The error message I've encountered is; ERROR:root:Error generating model response: exception: access violation writing 0x0000000000005140, which seems to be indicative of an access violation...
1
by: PapaRatzi | last post by:
Hello, I am teaching myself MS Access forms design and Visual Basic. I've created a table to capture a list of Top 30 singles and forms to capture new entries. The final step is a form (unbound)...
1
by: CloudSolutions | last post by:
Introduction: For many beginners and individual users, requiring a credit card and email registration may pose a barrier when starting to use cloud servers. However, some cloud server providers now...
0
by: af34tf | last post by:
Hi Guys, I have a domain whose name is BytesLimited.com, and I want to sell it. Does anyone know about platforms that allow me to list my domain in auction for free. Thank you

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.