令人惊讶的是之前没有人提出这个问题，但是对于需要数据库的实用主义者来说:MongoDB对Geo查询提供了出色的支持，包括这个查询。

你需要的是:

db.neighborhoods。findOne({geometry: {$geoIntersects: {$geometry: { type: "Point"，坐标:["经度"，"纬度"]}}} })

communities是存储一个或多个标准GeoJson格式多边形的集合。如果查询返回null，则表示不相交，否则为。

这里有详细的记录: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/

在330个不规则多边形网格中，超过6000个点分类的性能不到一分钟，没有任何优化，包括用各自的多边形更新文档的时间。

这个问题的大多数答案并没有很好地处理所有的极端情况。以下是一些微妙的极端情况: 这是一个javascript版本，所有角落的情况都得到了很好的处理。

/** Get relationship between a point and a polygon using ray-casting algorithm
 * @param {{x:number, y:number}} P: point to check
 * @param {{x:number, y:number}[]} polygon: the polygon
 * @returns -1: outside, 0: on edge, 1: inside
 */
function relationPP(P, polygon) {
    const between = (p, a, b) => p >= a && p <= b || p <= a && p >= b
    let inside = false
    for (let i = polygon.length-1, j = 0; j < polygon.length; i = j, j++) {
        const A = polygon[i]
        const B = polygon[j]
        // corner cases
        if (P.x == A.x && P.y == A.y || P.x == B.x && P.y == B.y) return 0
        if (A.y == B.y && P.y == A.y && between(P.x, A.x, B.x)) return 0

        if (between(P.y, A.y, B.y)) { // if P inside the vertical range
            // filter out "ray pass vertex" problem by treating the line a little lower
            if (P.y == A.y && B.y >= A.y || P.y == B.y && A.y >= B.y) continue
            // calc cross product `PA X PB`, P lays on left side of AB if c > 0 
            const c = (A.x - P.x) * (B.y - P.y) - (B.x - P.x) * (A.y - P.y)
            if (c == 0) return 0
            if ((A.y < B.y) == (c > 0)) inside = !inside
        }
    }

    return inside? 1 : -1
}

这只适用于凸形状，但是Minkowski Portal Refinement和GJK也是测试一个点是否在多边形中的很好的选择。您使用闵可夫斯基减法从多边形中减去点，然后运行这些算法来查看多边形是否包含原点。

另外，有趣的是，你可以用支持函数更隐式地描述你的形状，它以一个方向向量作为输入，并输出沿该向量的最远点。这可以让你描述任何凸形状..弯曲的，由多边形制成的，或混合的您还可以执行一些操作，将简单支持函数的结果组合起来，以生成更复杂的形状。

更多信息: http://xenocollide.snethen.com/mpr2d.html

此外，game programming gems 7讨论了如何在3d中做到这一点(:

没有什么比归纳定义问题更美好的了。为了完整起见，你在序言中有一个版本，它可能也澄清了光线投射背后的思想:

基于仿真的简化算法在http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

一些helper谓词:

exor(A,B):- \+A,B;A,\+B.
in_range(Coordinate,CA,CB) :- exor((CA>Coordinate),(CB>Coordinate)).

inside(false).
inside(_,[_|[]]).
inside(X:Y, [X1:Y1,X2:Y2|R]) :- in_range(Y,Y1,Y2), X > ( ((X2-X1)*(Y-Y1))/(Y2-Y1) +      X1),toggle_ray, inside(X:Y, [X2:Y2|R]); inside(X:Y, [X2:Y2|R]).

get_line(_,_,[]).
get_line([XA:YA,XB:YB],[X1:Y1,X2:Y2|R]):- [XA:YA,XB:YB]=[X1:Y1,X2:Y2]; get_line([XA:YA,XB:YB],[X2:Y2|R]).

给定两点a和B的直线(直线(a,B))方程为:

                    (YB-YA)
           Y - YA = ------- * (X - XA) 
                    (XB-YB) 

It is important that the direction of rotation for the line is setted to clock-wise for boundaries and anti-clock-wise for holes. We are going to check whether the point (X,Y), i.e the tested point is at the left half-plane of our line (it is a matter of taste, it could also be the right side, but also the direction of boundaries lines has to be changed in that case), this is to project the ray from the point to the right (or left) and acknowledge the intersection with the line. We have chosen to project the ray in the horizontal direction (again it is a matter of taste, it could also be done in vertical with similar restrictions), so we have:

               (XB-XA)
           X < ------- * (Y - YA) + XA
               (YB-YA) 

Now we need to know if the point is at the left (or right) side of the line segment only, not the entire plane, so we need to restrict the search only to this segment, but this is easy since to be inside the segment only one point in the line can be higher than Y in the vertical axis. As this is a stronger restriction it needs to be the first to check, so we take first only those lines meeting this requirement and then check its possition. By the Jordan Curve theorem any ray projected to a polygon must intersect at an even number of lines. So we are done, we will throw the ray to the right and then everytime it intersects a line, toggle its state. However in our implementation we are goint to check the lenght of the bag of solutions meeting the given restrictions and decide the innership upon it. for each line in the polygon this have to be done.

is_left_half_plane(_,[],[],_).
is_left_half_plane(X:Y,[XA:YA,XB:YB], [[X1:Y1,X2:Y2]|R], Test) :- [XA:YA, XB:YB] =  [X1:Y1, X2:Y2], call(Test, X , (((XB - XA) * (Y - YA)) / (YB - YA) + XA)); 
                                                        is_left_half_plane(X:Y, [XA:YA, XB:YB], R, Test).

in_y_range_at_poly(Y,[XA:YA,XB:YB],Polygon) :- get_line([XA:YA,XB:YB],Polygon),  in_range(Y,YA,YB).
all_in_range(Coordinate,Polygon,Lines) :- aggregate(bag(Line),    in_y_range_at_poly(Coordinate,Line,Polygon), Lines).

traverses_ray(X:Y, Lines, Count) :- aggregate(bag(Line), is_left_half_plane(X:Y, Line, Lines, <), IntersectingLines), length(IntersectingLines, Count).

% This is the entry point predicate
inside_poly(X:Y,Polygon,Answer) :- all_in_range(Y,Polygon,Lines), traverses_ray(X:Y, Lines, Count), (1 is mod(Count,2)->Answer=inside;Answer=outside).

下面是Rust版本的@nirg答案(Philipp Lenssen javascript版本) 我给出这个答案是因为我从这个网站得到了很多帮助，我翻译javascript版本rust作为一个练习，希望可以帮助一些人，最后一个原因是，在我的工作中，我会把这段代码翻译成一个wasm，以提高我的画布的性能，这是一个开始。我的英语很差……,请原谅我 ｀

pub struct Point {
    x: f32,
    y: f32,
}
pub fn point_is_in_poly(pt: Point, polygon: &Vec<Point>) -> bool {
    let mut is_inside = false;

    let max_x = polygon.iter().map(|pt| pt.x).reduce(f32::max).unwrap();
    let min_x = polygon.iter().map(|pt| pt.x).reduce(f32::min).unwrap();
    let max_y = polygon.iter().map(|pt| pt.y).reduce(f32::max).unwrap();
    let min_y = polygon.iter().map(|pt| pt.y).reduce(f32::min).unwrap();

    if pt.x < min_x || pt.x > max_x || pt.y < min_y || pt.y > max_y {
        return is_inside;
    }

    let len = polygon.len();
    let mut j = len - 1;

    for i in 0..len {
        let y_i_value = polygon[i].y > pt.y;
        let y_j_value = polygon[j].y > pt.y;
        let last_check = (polygon[j].x - polygon[i].x) * (pt.y - polygon[i].y)
            / (polygon[j].y - polygon[i].y)
            + polygon[i].x;
        if y_i_value != y_j_value && pt.x < last_check {
            is_inside = !is_inside;
        }
        j = i;
    }
    is_inside
}


let pt = Point {
    x: 1266.753,
    y: 97.655,
};
let polygon = vec![
    Point {
        x: 725.278,
        y: 203.586,
    },
    Point {
        x: 486.831,
        y: 441.931,
    },
    Point {
        x: 905.77,
        y: 445.241,
    },
    Point {
        x: 1026.649,
        y: 201.931,
    },
];
let pt1 = Point {
    x: 725.278,
    y: 203.586,
};
let pt2 = Point {
    x: 872.652,
    y: 321.103,
};
println!("{}", point_is_in_poly(pt, &polygon));// false
println!("{}", point_is_in_poly(pt1, &polygon)); // true
println!("{}", point_is_in_poly(pt2, &polygon));// true

`

VBA版本:

注意:请记住，如果你的多边形是地图中的一个区域，纬度/经度是Y/X值，而不是X/Y(纬度= Y，经度= X)，因为从我的理解来看，这是历史含义，因为经度不是一个测量值。

类模块:CPoint

Private pXValue As Double
Private pYValue As Double

'''''X Value Property'''''

Public Property Get X() As Double
    X = pXValue
End Property

Public Property Let X(Value As Double)
    pXValue = Value
End Property

'''''Y Value Property'''''

Public Property Get Y() As Double
    Y = pYValue
End Property

Public Property Let Y(Value As Double)
    pYValue = Value
End Property

模块:

Public Function isPointInPolygon(p As CPoint, polygon() As CPoint) As Boolean

    Dim i As Integer
    Dim j As Integer
    Dim q As Object
    Dim minX As Double
    Dim maxX As Double
    Dim minY As Double
    Dim maxY As Double
    minX = polygon(0).X
    maxX = polygon(0).X
    minY = polygon(0).Y
    maxY = polygon(0).Y

    For i = 1 To UBound(polygon)
        Set q = polygon(i)
        minX = vbMin(q.X, minX)
        maxX = vbMax(q.X, maxX)
        minY = vbMin(q.Y, minY)
        maxY = vbMax(q.Y, maxY)
    Next i

    If p.X < minX Or p.X > maxX Or p.Y < minY Or p.Y > maxY Then
        isPointInPolygon = False
        Exit Function
    End If


    ' SOURCE: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

    isPointInPolygon = False
    i = 0
    j = UBound(polygon)

    Do While i < UBound(polygon) + 1
        If (polygon(i).Y > p.Y) Then
            If (polygon(j).Y < p.Y) Then
                If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
                    isPointInPolygon = True
                    Exit Function
                End If
            End If
        ElseIf (polygon(i).Y < p.Y) Then
            If (polygon(j).Y > p.Y) Then
                If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
                    isPointInPolygon = True
                    Exit Function
                End If
            End If
        End If
        j = i
        i = i + 1
    Loop   
End Function

Function vbMax(n1, n2) As Double
    vbMax = IIf(n1 > n2, n1, n2)
End Function

Function vbMin(n1, n2) As Double
    vbMin = IIf(n1 > n2, n2, n1)
End Function


Sub TestPointInPolygon()

    Dim i As Integer
    Dim InPolygon As Boolean

'   MARKER Object
    Dim p As CPoint
    Set p = New CPoint
    p.X = <ENTER X VALUE HERE>
    p.Y = <ENTER Y VALUE HERE>

'   POLYGON OBJECT
    Dim polygon() As CPoint
    ReDim polygon(<ENTER VALUE HERE>) 'Amount of vertices in polygon - 1
    For i = 0 To <ENTER VALUE HERE> 'Same value as above
       Set polygon(i) = New CPoint
       polygon(i).X = <ASSIGN X VALUE HERE> 'Source a list of values that can be looped through
       polgyon(i).Y = <ASSIGN Y VALUE HERE> 'Source a list of values that can be looped through
    Next i

    InPolygon = isPointInPolygon(p, polygon)
    MsgBox InPolygon

End Sub