真的很喜欢Nirg发布的解决方案，由bobobobo编辑。我只是让它javascript友好，更容易读懂我的使用:

function insidePoly(poly, pointx, pointy) {
    var i, j;
    var inside = false;
    for (i = 0, j = poly.length - 1; i < poly.length; j = i++) {
        if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside;
    }
    return inside;
}

以下是M. Katz基于Nirg方法的答案的JavaScript变体:

function pointIsInPoly(p, polygon) {
    var isInside = false;
    var minX = polygon[0].x, maxX = polygon[0].x;
    var minY = polygon[0].y, maxY = polygon[0].y;
    for (var n = 1; n < polygon.length; n++) {
        var q = polygon[n];
        minX = Math.min(q.x, minX);
        maxX = Math.max(q.x, maxX);
        minY = Math.min(q.y, minY);
        maxY = Math.max(q.y, maxY);
    }

    if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
        return false;
    }

    var i = 0, j = polygon.length - 1;
    for (i, j; i < polygon.length; j = i++) {
        if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
                p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
            isInside = !isInside;
        }
    }

    return isInside;
}

下面是golang版本的@nirg答案(灵感来自于@@m-katz的c#代码)

func isPointInPolygon(polygon []point, testp point) bool {
    minX := polygon[0].X
    maxX := polygon[0].X
    minY := polygon[0].Y
    maxY := polygon[0].Y

    for _, p := range polygon {
        minX = min(p.X, minX)
        maxX = max(p.X, maxX)
        minY = min(p.Y, minY)
        maxY = max(p.Y, maxY)
    }

    if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
        return false
    }

    inside := false
    j := len(polygon) - 1
    for i := 0; i < len(polygon); i++ {
        if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
            inside = !inside
        }
        j = i
    }

    return inside
}

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

你需要的是:

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

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

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

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

Like David Segonds' answer suggests I use an approach of angle summation derived from my concave polygon drawing algorithm. It relies of adding up the approximate angles of subtriangles around the point to obtain a weight. A weight around 1.0 means the point is inside the triangle, a weight around 0.0 means outside, a weight around -1.0 is what happens when inside the polygon but in reverse order (like with one of the halves of a bowtie-shaped tetragon) and a weight of NAN if exactly on an edge. The reason it's not slow is that angles don't need to be estimated accurately at all. Holes can be handled by treating them as separate polygons and subtracting the weights.

typedef struct { double x, y; } xy_t;

xy_t sub_xy(xy_t a, xy_t b)
{
    a.x -= b.x;
    a.y -= b.y;
    return a;
}

double calc_sharp_subtriangle_pixel_weight(xy_t p0, xy_t p1)
{
    xy_t rot, r0, r1;
    double weight;

    // Rotate points (unnormalised)
    rot = sub_xy(p1, p0);
    r0.x = rot.x*p0.y - rot.y*p0.x;
    r0.y = rot.x*p0.x + rot.y*p0.y;
    r1.y = rot.x*p1.x + rot.y*p1.y;

    // Calc weight
    weight = subtriangle_angle_approx(r1.y, r0.x) - subtriangle_angle_approx(r0.y, r0.x);

    return weight;
}

double calc_sharp_polygon_pixel_weight(xy_t p, xy_t *corner, int corner_count)
{
    int i;
    xy_t p0, p1;
    double weight = 0.;

    p0 = sub_xy(corner[corner_count-1], p);
    for (i=0; i < corner_count; i++)
    {
        // Transform corner coordinates
        p1 = sub_xy(corner[i], p);

        // Calculate weight for each subtriangle
        weight += calc_sharp_subtriangle_pixel_weight(p0, p1);
        p0 = p1;
    }

    return weight;
}

因此，对于多边形的每一段，都形成一个子三角形，并计算点，然后旋转每个子三角形以计算其近似角度并添加到权重。

调用subtriangle_angle_approx(y, x)可以替换为atan2(y, x) / (2.*pi)，但是一个非常粗略的近似值就足够精确了:

double subtriangle_angle_approx(double y, double x)
{
    double angle, d;
    int obtuse;

    if (x == 0.)
        return NAN;

    obtuse = fabs(y) > fabs(x);
    if (obtuse)
        swap_double(&y, &x);

    // Core of the approximation, a very loosely approximate atan(y/x) / (2.*pi) over ]-1 , 1[
    d = y / x;
    angle = 0.13185 * d;

    if (obtuse)
        angle = sign(d)*0.25 - angle;

    return angle;
}

bobobobo引用的Eric Haines的文章真的很棒。特别有趣的是比较算法性能的表格;角度求和法和其他方法比起来真的很差。同样有趣的是，使用查找网格将多边形进一步细分为“in”和“out”扇区的优化可以使测试非常快，即使是在> 1000条边的多边形上。

不管怎样，现在还为时过早，但我的投票倾向于“交叉”方法，我认为这几乎就是Mecki所描述的。然而，我发现大卫·伯克(David Bourke)对它进行了最简洁的描述和编纂。我喜欢它不需要真正的三角函数，它适用于凸和凹，而且随着边数的增加，它的表现也相当不错。

顺便说一下，这是Eric Haines文章中的一个性能表，在随机多边形上进行测试。

                       number of edges per polygon
                         3       4      10      100    1000
MacMartin               2.9     3.2     5.9     50.6    485
Crossings               3.1     3.4     6.8     60.0    624
Triangle Fan+edge sort  1.1     1.8     6.5     77.6    787
Triangle Fan            1.2     2.1     7.3     85.4    865
Barycentric             2.1     3.8    13.8    160.7   1665
Angle Summation        56.2    70.4   153.6   1403.8  14693

Grid (100x100)          1.5     1.5     1.6      2.1      9.8
Grid (20x20)            1.7     1.7     1.9      5.7     42.2
Bins (100)              1.8     1.9     2.7     15.1    117
Bins (20)               2.1     2.2     3.7     26.3    278