Wednesday, December 28, 2011

Linking a point to the nearest line using OGR with Oracle Spatial

An Oracle table contains point data. Another contains line data. Both were imported from some CAD application that did not deem it important to consign the relationships between both. That was left to the programmer as an exercise.

Now, of course, using PL/SQL and Oracle spatial functions or operators is not an option. Why not? Because Oracle is not a good spatial RDBMS and PL/SQL is not a good scripting language. Maybe you could do it if you were willing to unlearn how to do things in a logical, straightforward GIS way and were willing to modify your data in a way to accommodate the idiosyncrasies of a crappy RDBMS (assuming you have admin powers on the Oracle server), but very few of you actually would. Unless, of course, you don't know better - in which case it's time to get wiser and learn something new. OGR and Python, for example, a magical pair, thick as thieves.

The solution is so simple, I'm ashamed to give it:
  def getClosestLine(point_geom):   
   mind = {}   
   for i in line_geoms:   
    d = point_geom.Distance(i)   
    mind[d] = i   
   return mind[min(mind.keys()]   

Nice solution but a bit slow. A quicker alternative, although less precise is to use incremental point buffering and selecting the first line it intersects, like this:
 def getClosestLine(point_geom):  
   found = False  
   buf_size = 0.0  
   step = 0.01  
   closest = None  
   while not found:  
     buf = point_geom.Buffer(buf_size)  
     intersecting = [ x for x in line_geoms if x.Intersect(buf) ]  
     if len(intersecting) > 0:  
       found = True  
       closest = intersecting[0] #can be multiple, not 100% precise  
       break  
     else:  
       buf_size += step  
     if buf_size > 10.0: break  
   return closest  

It's quicker than the previous (although still very slow), but it's not quite precise, as there will always be that margin for error the size of step. The finer the step you choose, the slower it'll get. But if you're only looking for an approximation, I guess it's OK.

Since I had a rough idea of the distance in which the lines were likely to occur (by testing), I could experience with a blend of both solutions. I created a rough buffer (1m) around my point and filtered the lines intersecting with it. I then calculated the minimum distance for only those lines to select my winner.
  def getClosestLine(point_geom):   
   mind = {}   
   buf_size = 1  
   buf = geom.Buffer(buf_size)  
   intersecting = [ x for x in line_geoms if x.Intersect(buf) ]  
   for i in intersecting:   
    d = point_geom.Distance(i)   
    mind[d] = i  
   return mind[min(mind.keys()]   

This was largely efficient for the tens of thousands of points and lines I had to relate. It kept me away from the coffee machine during testing and got me the bad eye from the Oracle boys too stubborn to learn something new and less expensive. Just kidding, the bad eye wasn't directed towards me - it's just the glum stare you get when working too long with bad software! Ha ha

No comments:

Post a Comment