[FWTools] GEOS C API does not include DistanceOp?

[Bryan Keith]

Frank Warmerdam wrote:
> Bryan Keith wrote:
>> Hello all,
>>
>> I'm using FWTools1.0.0 on windows.  I'm trying to get the distance from 
>> a point to a polygon (or line), but I get this error using the python 
>> wrapper:
>>
>> GEOS C API does not include DistanceOp
>>
>> Is there a good workaround?  Anyone have a Python algorithm to get the 
>> distance between ogr geometries?
>>
>> I also have FWTools-linux-1.0.0a7 installed on a linux machine, but I 
>> get the same GEOS error.
> 
> Bryan,
> 
> If you file a bug against GDAL/OGR on this, I'll look at correcting the
> problem.  The root of the problem was that the GEOS Distance operation was
> not originally exposed in the GEOS C API, so when I shifted to that it
> was lost.  I gather it is now available in the latest GEOS releases, so
> the problem should be resolvable.

Frank,

Thanks.  I added a bug and assigned it to you:

http://bugzilla.remotesensing.org/show_bug.cgi?id=1227

> 
> I'm afraid I don't have a good workaround suggestion in the meantime.

In the meantime I wrote some Python code that seems to be working for my 
specific situation:  one of my geometries is always a simple point and 
I'm assuming everything is on the same z plane.  I'm attaching the code 
to this e-mail if anyone's interested.  Comments are appreciated.

Bryan

> 
> Best regards,
-------------- next part --------------
import ogr
def ReturnPTLNDistance(mPT, mLN, dDist = "None"):
  #mPT is an ogr point
  #mLN is an ogr polygon or polyline
  #dDist is the current calculated minimum distance, used when this function
  #calls itself; leave dDist out when calling from another function (when no
  #minimum distance is yet determined)
  px = mPT.GetX()
  py = mPT.GetY()
  if (mLN.GetGeometryCount() == 0):
    for j in range(mLN.GetPointCount() - 1):
      X1 = mLN.GetX(j)
      Y1 = mLN.GetY(j)
      X2 = mLN.GetX(j + 1)
      Y2 = mLN.GetY(j + 1)
      d = DistancePointLine(px, py, X1, Y1, X2, Y2)
      if (dDist == "None"):
        dDist = d
      else:
        dDist = min(dDist, d)
  else:
    for i in range(mLN.GetGeometryCount()):
      myGeom = mLN.GetGeometryRef(i)
      dDist = ReturnPTLNDistance(mPT, myGeom, dDist)
  return dDist
def ReturnDistance(x1, x2, y1, y2, z1 = 0, z2 = 0):
  return ((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2) ** 0.5
def DistancePointLine(px, py, x1, y1, x2, y2):
  #px,py is the point to test.
  #x1,y1,x2,y2 is the line to check distance.
  #Returns distance from the line, or if the intersecting point on the line
  #nearest the point tested is outside the endpoints of the line, the distance to
  #the nearest endpoint.
  #code taken from VB code from here:
  #http://astronomy.swin.edu.au/~pbourke/geometry/pointline/source.vba
  LineMag = ReturnDistance(x1, x2, y1, y2)
  if (LineMag < 0.00000001):
    return ReturnDistance(x1, px, y1, py)
  u = ((px - x1) * (x2 - x1)) + ((py - y1) * (y2 - y1))
  u = u / (LineMag ** 2)
  if ((u < 0.00001) or (u > 1)):
    #closest point does not fall within the line segment, take the shorter
    #distance to an endpoint
    ix = ReturnDistance(px, x1, py, y1)
    iy = ReturnDistance(px, x2, py, y2)
    if (ix > iy):
      return iy
    else:
      return ix
  else:
    #Intersecting point is on the line, use the formula
    #see http://astronomy.swin.edu.au/~pbourke/geometry/pointline/
    ix = x1 + u * (x2 - x1)
    iy = y1 + u * (y2 - y1)
    return ReturnDistance(px, ix, py, iy)