Wednesday, December 28, 2011

Getting a point on a line at a certain distance using OGR

Another day, another adventure. I arrived at work and was assigned the task to look for a point at a certain distance on a line. The point didn't exist yet, I had to create it and return it. My boss, grumpy and suffering from a hangover, was threatening (like the whole week before that) to get my ass fired if I didn't prove myself more productive and resourceful. So I grabbed myself a mug of coffee and set myself on the search of that elusive point, having yet to pay the month's rent and put the food on the table for the wife and cats.

The distance was relative: it was to be the middle of the segment of the line where the middle of the line occurred. Was that clear? A line can have several segments. The segment that interests us is the one on which the line hits half of its total distance. We want the point at the middle of that segment. We want it real bad!



A quick look at the OGR API tells us there's a Value method for OGRLineString. But that's not exactly what we want, and there's no Python binding anyway. So, God, what am I gonna do?

"Don't panic." a calm voice told me. "You have good friends and they've dragged you out of a bigger mess than this one before."

So I calmed down and started coding, and very soon I was feeling better. Here's is the piece of code that flowed from my fingertips:
 for line_feat in layer:  
   geom = line_feat.geometry()  
   num_p = geom.GetPointCount()  
   current_distance = 0.0  
   half_distance = geom.Length() / 2  
   for i in xrange(num_p-1):  
     if i+1 < num_p:      
       line = ogr.Geometry(ogr.wkbLineString)  
       line.AddPoint(geom.GetPoint(i)[0],geom.GetPoint(i)[1])  
       line.AddPoint(geom.GetPoint(i+1)[0],geom.GetPoint(i+1)[1])  
       segment_len = line.Length()        
       if (current_distance + segment_len) > half_distance:  
         x_0 = line.GetPoint(0)[0]  
         y_0 = line.GetPoint(0)[1]  
         x_1 = line.GetPoint(1)[0]  
         y_1 = line.GetPoint(1)[1]  
         d_x = (x_1 - x_0) / 2  
         d_y = (y_1 - y_0) / 2  
         x_n = x_1 - d_x  
         y_n = y_1 - d_y          
         createPoint(x_n,y_n)  
         break  
       current_distance = current_distance + segment_len  

In short:
  • Get the half-length of the total line feature
  • Obtain the points from which you reconstitute each segment composing the line
  • You calculate the length of these segments and sum them up, when you've exceeded the half-length of the total line, you know you're holding your precious segment
  • You calculate the point at the middle of that segment - that's your point!

No comments:

Post a Comment