distance from point to polygon (when inside)

Question

I'm trying to find an efficient way to compute the distance from a point to the nearest edge of a polygon in python. I thought shapely would be perfect for this, but it only computes the distance when the point is outside the polygon. I thought about opencv, which has a built in function, but it requires integer coordinates (I can't discretize my data). The polygon is not guaranteed to be convex.

Answer 1

Use the .boundary property of the Polygon to return a LineString of the LinearRings, then use the .distance property. You can also use the .exterior property if you want only the exterior ring, and not any of the interior rings (if they exist). E.g.

from shapely import wkt
poly = wkt.loads('POLYGON((30 10, 40 40, 20 40, 10 20, 30 10))')
pt = wkt.loads('POINT(20 20)')

# The point is in the polygon, so the distance will always be 0.0
poly.distance(pt)  # 0.0

# The distance is calculated to the nearest boundary
poly.boundary.distance(pt)  # 4.47213595499958

# or exterior ring
poly.exterior.distance(pt)  # 4.47213595499958

The result is the same for .boundary and .exterior with this example, but could change for other examples.

Answer 2

I have some code to get line and polygon intersection (given by edges) in cython. If You dont want to wait here it is. Im going to sleep now, so its only raw code of ship hull boundary and waterline intersection. Have to implement whats needed on Your own for today :P...

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.interpolate as sci_itp
import scipy.integrate as sci_itg
cimport numpy as np
cimport cython

from CROSS_SECTIONS_universal import *

DTYPE = np.float64
ctypedef np.float64_t DTYPE_t

ITYPE = np.int
ctypedef np.int_t ITYPE_t

ITYPEP = np.intp
ctypedef np.intp_t ITYPEP_t

ITYPE1 = np.int64
ctypedef np.int64_t ITYPE1_t

BTYPE = np.uint8
ctypedef np.uint8_t BTYPE_t

ctypedef Py_ssize_t shape_t

BOTYPE = np.bool_

cdef double[:,:] waterLine = np.array([[0,0],[0,0]],dtype=DTYPE)
cdef double[:,:] crossLine = np.array([[0,0],[0,0]],dtype=DTYPE)
cdef double[:,:] line
cdef double[:] sectionLinePoint = np.full((2),np.nan,dtype=DTYPE)
cdef double[:,:] sectionLengthPOINTS = np.full((2,2),np.nan,dtype=DTYPE)
cdef double[:,:] sectionBreadthPOINTS = np.full((2,2),np.nan,dtype=DTYPE)
cdef double sectionLength,sectionBreadth


cdef point_side_of_line(double ax,double ay,double bx,double by,double cx,double cy):
    """ Returns a position of the point c relative to the line going through a and b
        Points a, b are expected to be different
    """
    cdef double s10x,s10y,s20x,s20y,denom
    cdef bint denomIsPositive,denomIsNegative
    cdef int side

    s10x,s10y=bx-ax,by-ay
    s20x,s20y=cx-ax,cy-ay

    denom = s20y*s10x - s10y*s20x
    denomIsPositive,denomIsNegative=denom>0,denom<0

    side=1 if denomIsPositive else (-1 if denomIsNegative else 0)
    return side

cdef is_point_in_closed_segment(double ax,double ay,double bx,double by,double cx,double cy):
    """ Returns True if c is inside closed segment, False otherwise.
        a, b, c are expected to be collinear
    """
    cdef bint pointIn

    if ax < bx:
        pointIn=ax <= cx and cx <= bx
        return pointIn
    if bx < ax:
        pointIn=bx <= cx and cx <= ax
        return pointIn
    if ay < by:
        pointIn=ay <= cy and cy <= by
        return pointIn
    if by < ay:
        pointIn=by <= cy and cy <= ay
        return pointIn

    pointIn=ax == cx and ay == cy
    return pointIn

cdef closed_segment_intersect(double ax,double ay,double bx,double by,double cx,double cy,double dx,double dy):
    """ Verifies if closed segments a, b, c, d do intersect.
    """
    cdef int s1,s2

    if (ax is bx) and (ay is by):
        assert (ax is cx and ay is cy) or (ax is dx and ay is dy)
    if (cx is dx) and (cy is dy):
        assert (cx is ax and cy is ay) or (cx is bx and cy is by)

    s1,s2 = point_side_of_line(ax,ay,bx,by,cx,cy),point_side_of_line(ax,ay,bx,by,dx,dy)

    # All points are collinear
    if s1 is 0 and s2 is 0:
        assert \
            is_point_in_closed_segment(ax,ay,bx,by,cx,cy) or \
            is_point_in_closed_segment(ax,ay,bx,by,dx,dy) or \
            is_point_in_closed_segment(cx,cy,dx,dy,ax,ay) or \
            is_point_in_closed_segment(cx,cy,dx,dy,bx,by)


    # No touching and on the same side
    if s1 and s1 is s2:
        assert False

    s1,s2 = point_side_of_line(cx,cy,dx,dy,ax,ay),point_side_of_line(cx,cy,dx,dy,bx,by)

    # No touching and on the same side
    if s1 and s1 is s2:
        assert False

cdef find_intersection(double[:,:] L1,double[:,:] L2, double[:] intersectionPoint) :
    cdef double ax,ay,bx,by,cx,cy,dx,dy
    cdef double s10x,s10y,s32x,s32y,s02x,s02y
    cdef double tNumer,t,denom


    ax,ay=L1[0]
    bx,by=L1[1]
    cx,cy=L2[0]
    dx,dy=L2[1]

    closed_segment_intersect(ax,ay,bx,by,cx,cy,dx,dy)

    s10x,s10y=bx-ax,by-ay
    s02x,s02y=ax-cx,ay-cy
    s32x,s32y=dx-cx,dy-cy

    denom = s10x * s32y - s32x * s10y
    tNumer = s32x * s02y - s32y * s02x

    t = tNumer / denom
    intersectionPoint[0]=ax+t*s10x
    intersectionPoint[1]=ay+t*s10y

cdef section_length_breadth(double[:,:] planeSectionCLUSTER,double x0,double x1, double y0, double y1):

    cdef int counterL=0,counterB=0

    x01=(x0+x1)/2.

    crossLine[0,0]=x01
    crossLine[1,0]=x01
    crossLine[0,1]=y0-0.1
    crossLine[1,1]=y1+0.1
    waterLine[0,0]=x0-0.1
    waterLine[1,0]=x1+0.1


    for i in range(1,len(planeSectionCLUSTER)):
        line=planeSectionCLUSTER[i-1:i+1]
        plt.plot(line[:,0],line[:,1],'c-',lw=3)
        try:
            try:
                find_intersection(line,waterLine,sectionLinePoint)
                assert sectionLinePoint[0] not in sectionLengthPOINTS[:,0]
                sectionLengthPOINTS[counterL]=sectionLinePoint
                counterL+=1
            except AssertionError:          
                find_intersection(line,crossLine,sectionLinePoint)
                assert sectionLinePoint[1] not in sectionBreadthPOINTS[:,1]
                sectionBreadthPOINTS[counterB]=sectionLinePoint
                counterB+=1
        except AssertionError:
            pass    

    print '#'
    print [[j for j in i] for i in sectionLengthPOINTS]
    print [[j for j in i] for i in sectionBreadthPOINTS]

    plt.plot(waterLine[:,0],waterLine[:,1],'b--',lw=3)
    plt.plot(sectionLengthPOINTS[:,0],sectionLengthPOINTS[:,1],'ko-',lw=3)
    plt.plot(crossLine[:,0],crossLine[:,1],'b--',lw=3)
    plt.plot(sectionBreadthPOINTS[:,0],sectionBreadthPOINTS[:,1],'ro-',lw=3)
    return 0,0