点到多边形的距离（在内部时）

Question

我正在尝试寻找一种有效的方法来计算 python 中一个点到多边形最近边的距离。我认为 shapely 非常适合这个，但它只计算点在多边形之外时的距离。我想到了 opencv，它有一个内置函数，但它需要整数坐标（我不能离散化我的数据）。不保证多边形是凸的。

Answer 1

用Polygon的.boundary属性到returnLinearRings的LineString，再用.distance属性。如果只需要外环而不需要任何内环（如果存在），也可以使用 .exterior 属性。例如

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

此示例的 .boundary 和 .exterior 的结果相同，但其他示例可能会有所不同。

Answer 2

我有一些代码可以在 cython 中获取线和多边形的交点（由边给出）。如果您不想在这里等待，那就是。我现在要睡觉了，所以它只是船体边界和水线交点的原始代码。必须自己实施今天需要的东西 :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

点到多边形的距离（在内部时）

distance from point to polygon (when inside)

python

geometry

numpy

point-in-polygon

shapely