确定球体上的经纬度矩形和圆是否重叠

发布于 2024-07-10 16:07:41 字数 218 浏览 7 评论 0原文

假设我有以下内容:

  • 由最小和最大纬度和经度定义的区域(通常是“经纬度矩形”,尽管除了某些投影之外它实际上不是矩形)。
  • 由中心纬度/经度和半径定义的圆

如何确定:

  1. 两个形状是否重叠?
  2. 圆是否完全包含在矩形内?

我正在寻找一个完整的公式/算法,而不是数学课本身。

Suppose I have the following:

  • A region defined by minimum and maximum latitude and longitude (commonly a 'lat-long rect', though it's not actually rectangular except in certain projections).
  • A circle, defined by a center lat/long and a radius

How can I determine:

  1. Whether the two shapes overlap?
  2. Whether the circle is entirely contained within the rect?

I'm looking for a complete formula/algorithm, rather than a lesson in the math, per-se.

如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。

评论(7

阳光下慵懒的猫 2024-07-17 16:07:41

警告:如果圆/“矩形”跨越球体的大部分,例如:

“矩形”:最小长 = -90deg,最大长 = +90deg,最小纬度 = +70deg ,最大纬度 = +80 度

圆:中心 = 纬度 = +85 度,长 = +160 度,半径 = 20 度(例如,如果 A 点在圆上,C 点是圆的中心,O 点是球体的中心,则角度AOC = 40 度)。

它们相交,但数学上可能有几种情况来检查相交/包含。 以下点位于上述圆上:P1=(+65 度纬度,+160 度长),P2=(+75 度纬度,-20 度长)。 P1 在“矩形”外部,P2 在“矩形”内部,因此圆/“矩形”至少有 2 个点相交。

好的,这是我对解决方案的概述:


设 C = 半径为 R 的圆心(如上所示表示为球面角)。 C 具有纬度 LATC 和经度 LONGC。 由于“矩形”这个词在这里有点误导(恒定纬度的线不是大圆的线段),我将使用术语“边界框”。

  • 函数 InsideCircle(P) 返回 +1、0 或 -1: 如果点 P 在圆内,则返回 +1;如果点 P 在圆上,则返回 0 ,如果 P 点在圆外,则为 -1:计算 C 和任意点 P 之间的大圆距离 D(表示为球面角)将告诉您 P 是否在圆内: InsideCircle(P) =sign(RD) (正如 Sente 中的用户 @Die 提到的,在这个论坛的其他地方已经询问过大圆距离)

  • 定义PANG(x) = x 的主角 = MOD(x+180deg, 360deg)-180deg。 PANG(x) 始终在 -180deg 和 +180deg 之间(包含 -180deg 和 +180deg)(+180deg 应映射到 -180deg)。

  • 要定义边界框,您需要知道 4 个数字,但经度有一个小问题。 LAT1和LAT2表示边界纬度(假设LAT1<LAT2); 那里没有任何歧义。 LONG1 和 LONG2 表示经度间隔的边界经度,但这会很棘手,并且更容易将此间隔重写为中心和宽度,其中 LONGM = 该间隔的中心,LONGW = 宽度。 请注意,经度间隔始终有 2 种可能性。 您必须指定是哪种情况,无论是包括还是排除 180 度子午线,例如,从 -179 度到 +177 度的最短间隔为 LONGM = +179 度和 LONGW = 4 度,但从 -179 度到 + 的其他间隔177deg 的 LONGM = -1deg 且 LONGW = 356deg。 如果您天真地尝试与区间 [-179,177] 进行“常规”比较,您最终将使用更大的区间,而这可能不是您想要的。 顺便说一句,如果以下两个条件都成立,则具有纬度 LATP 和经度 LONGP 的点 P 位于边界框内:

    • LAT1 <= LATP 和 LATP <= LAT2(这部分是显而易见的)
    • abs(PANG(LONGP-LONGM)) <; 长/2

如果 PTEST = union(PCORNER,PLAT,PLONG) 中的以下任何点 P 如下所述,则圆与边界框相交,但 InsideCircle() 并不都返回相同的结果:

  • PCORNER = 边界框的 4 个角
  • 点 PLAT 位于边界框两侧(没有或有 2 个),它们与圆心具有相同的纬度,如果 LATC 位于 LAT1 和 LAT2 之间,其中这些点的纬度为 LATC,经度为 LONG1 和 LONG2。
  • 边界框两侧的 PLONG 点(要么没有,要么有 2 个或 4 个!),它们与圆的中心具有相同的经度。 这些点的经度 EITHER = LONGC OR 经度 PANG(LONGC-180)。 如果abs(PANG(LONGC-LONGM)) < LONGW/2 则 LONGC 是有效的经度。 如果abs(PANG(LONGC-180-LONGM)) < LONGW/2 则 PANG(LONGC-180) 是有效经度。 这些经度中的一个或两个或都不在边界框的经度区间内。 选择具有有效经度、纬度 LAT1 和 LAT2 的点 PLONG。

上面列出的这些点 PLAT 和 PLONG 是边界框上距离圆“最近”的点(如果角不是;我在引号中使用“最近”,从纬度/长距离的意义上来说,不是很大-圆距离),并涵盖圆的中心位于边界框边界的一侧但圆上的点“潜行穿过”边界框边界的情况。

如果 PTEST 中的所有点 P 返回 InsideCircle(P) == +1(全部在圆内),则该圆包含整个边界框。

如果 PTEST 中的所有点 P 返回 InsideCircle(P) == -1(全部在圆之外),则圆完全包含在边界框内。

否则,圆和边界框之间至少有一个交点。 请注意,这不会计算这些点的位置,尽管如果您在 PTEST 中取任意 2 个点 P1 和 P2,其中 InsideCircle(P1) = -InsideCircle(P2),那么您可以通过二等分找到交点(效率低下)。 (如果 InsideCircle(P) 返回 0,则有一个交点,尽管浮点数学中的等式通常不可信。)

可能有一种更有效的方法可以做到这一点,但上述方法应该可行。

warning: this can be tricky if the circles / "rectangles" span large portions of the sphere, e.g.:

"rectangle": min long = -90deg, max long = +90deg, min lat = +70deg, max lat = +80deg

circle: center = lat = +85deg, long = +160deg, radius = 20deg (e.g. if point A is on the circle and point C is the circle's center, and point O is the sphere's center, then angle AOC = 40deg).

These intersect but the math is likely to have several cases to check intersection/containment. The following points lie on the circle described above: P1=(+65deg lat,+160deg long), P2=(+75deg lat, -20deg long). P1 is outside the "rectangle" and P2 is inside the "rectangle" so the circle/"rectangle" intersect in at least 2 points.

OK, here's my shot at an outline of the solution:


Let C = circle center with radius R (expressed as a spherical angle as above). C has latitude LATC and longitude LONGC. Since the word "rectangle" is kind of misleading here (lines of constant latitude are not segments of great circles), I'll use the term "bounding box".

  • function InsideCircle(P) returns +1,0,or -1: +1 if point P is inside the circle, 0 if point P is on the circle, and -1 if point P is outside the circle: calculation of great-circle distance D (expressed as spherical angle) between C and any point P will tell you whether or not P is inside the circle: InsideCircle(P) = sign(R-D) (As user @Die in Sente mentioned, great circle distances have been asked on this forum elsewhere)

  • Define PANG(x) = the principal angle of x = MOD(x+180deg, 360deg)-180deg. PANG(x) is always between -180deg and +180deg, inclusive (+180deg should map to -180deg).

  • To define the bounding box, you need to know 4 numbers, but there's a slight issue with longitude. LAT1 and LAT2 represent the bounding latitudes (assuming LAT1 < LAT2); there's no ambiguity there. LONG1 and LONG2 represent the bounding longitudes of a longitude interval, but this gets tricky, and it's easier to rewrite this interval as a center and width, with LONGM = the center of that interval and LONGW = width. NOTE that there are always 2 possibilities for longitude intervals. You have to specify which of these cases it is, whether you are including or excluding the 180deg meridian, e.g. the shortest interval from -179deg to +177deg has LONGM = +179deg and LONGW = 4deg, but the other interval from -179deg to +177deg has LONGM = -1deg and LONGW = 356deg. If you naively try to do "regular" comparisons with the interval [-179,177] you will end up using the larger interval and that's probably not what you want. As an aside, point P, with latitude LATP and longitude LONGP, is inside the bounding box if both of the following are true:

    • LAT1 <= LATP and LATP <= LAT2 (that part is obvious)
    • abs(PANG(LONGP-LONGM)) < LONGW/2

The circle intersects the bounding box if ANY of the following points P in PTEST = union(PCORNER,PLAT,PLONG) as described below, do not all return the same result for InsideCircle():

  • PCORNER = the bounding box's 4 corners
  • the points PLAT on the bounding box's sides (there are either none or 2) which share the same latitude as the circle's center, if LATC is between LAT1 and LAT2, in which case these points have the latitude LATC and longitude LONG1 and LONG2.
  • the points PLONG on the bounding box's sides (there are either none or 2 or 4!) which share the same longitude as the circle's center. These points have EITHER longitude = LONGC OR longitude PANG(LONGC-180). If abs(PANG(LONGC-LONGM)) < LONGW/2 then LONGC is a valid longitude. If abs(PANG(LONGC-180-LONGM)) < LONGW/2 then PANG(LONGC-180) is a valid longitude. Either or both or none of these longitudes may be within the longitude interval of the bounding box. Choose points PLONG with valid longitudes, and latitudes LAT1 and LAT2.

These points PLAT and PLONG as listed above are the points on the bounding box that are "closest" to the circle (if the corners are not; I use "closest" in quotes, in the sense of lat/long distance and not great-circle distance), and cover the cases where the circle's center lies on one side of the bounding box's boundary but points on the circle "sneak across" the bounding box boundary.

If all points P in PTEST return InsideCircle(P) == +1 (all inside the circle) then the circle contains the bounding box in its entirety.

If all points P in PTEST return InsideCircle(P) == -1 (all outside the circle) then the circle is contained entirely within the bounding box.

Otherwise there is at least one point of intersection between the circle and the bounding box. Note that this does not calculate where those points are, although if you take any 2 points P1 and P2 in PTEST where InsideCircle(P1) = -InsideCircle(P2), then you could find a point of intersection (inefficiently) by bisection. (If InsideCircle(P) returns 0 then you have a point of intersection, though equality in floating-point math is generally not to be trusted.)

There's probably a more efficient way to do this but the above should work.

一直在等你来 2024-07-17 16:07:41

使用立体投影。 所有圆(特别是纬度、经度和您的圆)都映射到平面中的圆(或线)。 现在这只是一个关于平面几何中的圆和线的问题(更好的是,所有的经线都是经过 0 的线,所有的纬度都是围绕 0 的圆)

Use the Stereographic projection. All circles (specifically latitudes, longitudes and your circle) map to circles (or lines) in the plane. Now it's just a question about circles and lines in plane geometry (even better, all the longitues are lines through 0, and all the latitudes are circles around 0)

满地尘埃落定 2024-07-17 16:07:41
  • 是的,如果盒子的角包含圆心。
  • 是的,如果任何盒子角位于圆心半径内。
  • 是的,如果盒子包含圆心的经度,并且最接近圆心纬度的盒子纬度的经度交点在圆心的半径内。
  • 是的,如果盒子包含圆心的纬度,并且最短交点方位距圆心的半径距离处的点“超出”最近的盒子经度; 其中最短交点方位角是通过查找从圆心到纬度为零且经度“超出”最近的框经度的点的初始方位角来确定的。
  • 不,否则。

假设:

  • 您可以找到从 A 点到 B 点的最小航向的初始方位。
  • 您可以找到两点之间的距离。

第一次检查很简单。 第二次检查只需要找到四个距离。 第三次检查只需要找到从圆心到(最近框纬度,圆心经度)的距离。

第四个检查需要找到最接近圆心的边界框的经度线。 然后找到距离圆心最远的经线所在的大圆的中心。 求从圆心到大圆心的初始方位角。 找到该轴承上距圆心的点圆半径。 如果该点位于距离圆中心最近的经度线的另一侧,则圆和边界框在该侧相交。

我觉得这应该是有缺陷的,但是我一直没能发现。

我似乎无法解决的真正问题是找到完美包含圆的边界框(对于不包含极点的圆)。 最小/最大纬度的方位似乎是圆心和圆半径/(球体周长/4)纬度的函数。 在赤道附近,它下降到 pi/2(东)或 3*pi/2(西)。 当中心接近极点且半径接近球周长/4 时,方位角接近零(北)或 pi(南)。

  • Yes, if the box corners contain the circle-center.
  • Yes, if any of the box corners are within radius of circle-center.
  • Yes, if the box contains the longitude of circle-center and the longitude intersection of the box-latitude closest to circle-center-latitude is within radius of circle-center.
  • Yes, if the box contains the latitude of circle-center and the point at radius distance from circle-center on shortest-intersection-bearing is "beyond" the closest box-longitude; where shortest-intersection-bearing is determined by finding the initial bearing from circle-center to a point at latitude zero and a longitude that is pi/2 "beyond" the closest box-longitude.
  • No, otherwise.

Assumptions:

  • You can find the initial-bearing of a minimum course from point A to point B.
  • You can find the distance between two points.

The first check is trivial. The second check just requires finding the four distances. The third check just requires finding the distance from circle-center to (closest-box-latitude, circle-center-longitude).

The fourth check requires finding the longitude line of the bounding box that is closest to the circle-center. Then find the center of the great circle on which that longitude line rests that is furthest from circle-center. Find the initial-bearing from circle-center to the great-circle-center. Find the point circle-radius from circle-center on that bearing. If that point is on the other side of the closest-longitude-line from circle-center, then the circle and bounding box intersect on that side.

It seems to me that there should be a flaw in this, but I haven't been able to find it.

The real problem that I can't seem to solve is to find the bounding-box that perfectly contains the circle (for circles that don't contain a pole). The bearing to the latitude min/max appears to be a function of the latitude of circle-center and circle-radius/(sphere circumference/4). Near the equator, it falls to pi/2 (east) or 3*pi/2 (west). As the center approaches the pole and the radius approaches sphere-circumference/4, the bearing approach zero (north) or pi (south).

温暖的光 2024-07-17 16:07:41

这个怎么样?

找到连接矩形中心(点 Cr)和圆心的向量 v。 找到 v 与矩形相交的点 i。 如果<代码>||i-Cr|| + r > ||v|| 然后它们相交。

换句话说,矩形内部线段的长度加上圆形内部线段的长度应该大于(v,中心连接线段)的总长度。

找到点i应该是棘手的部分,特别是如果它落在经度边缘上,但你应该能够比我更快地找到一些东西。

编辑:此方法无法判断圆是否完全在矩形内。 为此,您需要找到从其中心到矩形所有四个边缘的距离。

编辑:上面是不正确的。正如 Federico Ramponi 所建议的那样,在某些情况下,即使在欧几里得几何中它也不起作用。 我会发布另一个答案。 请不接受这一点并随时投反对票。 我很快就会删除它。

How about this?

Find vector v that connects the center of the rectangle, point Cr, to the center of the circle. Find point i where v intersects the rectangle. If ||i-Cr|| + r > ||v|| then they intersect.

In other words, the length of the segment inside the rectangle plus the length of the segment inside the circle should be greater than the total length (of v, the center-connecting line segment).

Finding point i should be the tricky part, especially if it falls on a longitude edge, but you should be able to come up with something faster than I can.

Edit: This method can't tell if the circle is completely within the rectangle. You'd need to find the distance from its center to all four of the rectangle's edges for that.

Edit: The above is incorrect. There are some cases, as Federico Ramponi has suggested, where it does not work even in Euclidean geometry. I'll post another answer. Please unaccept this and feel free to vote down. I'll delete it shortly.

随遇而安 2024-07-17 16:07:41

这应该适用于地球上的任何点。 如果您想将其更改为不同大小的球体,只需将 kEarchRadiusKms 更改为您想要的球体半径即可。

该方法用于计算纬度点和经度点之间的距离。

我从这里得到了这个距离公式:
http://www.codeproject.com/csharp/distance Betweenlocations.asp

public static double Calc(double Lat1, double Long1, double Lat2, double Long2)
{
    double dDistance = Double.MinValue;
    double dLat1InRad = Lat1 * (Math.PI / 180.0);
    double dLong1InRad = Long1 * (Math.PI / 180.0);
    double dLat2InRad = Lat2 * (Math.PI / 180.0);
    double dLong2InRad = Long2 * (Math.PI / 180.0);

    double dLongitude = dLong2InRad - dLong1InRad;
    double dLatitude = dLat2InRad - dLat1InRad;

    // Intermediate result a.
    double a = Math.Pow(Math.Sin(dLatitude / 2.0), 2.0) +
               Math.Cos(dLat1InRad) * Math.Cos(dLat2InRad) *
               Math.Pow(Math.Sin(dLongitude / 2.0), 2.0);

    // Intermediate result c (great circle distance in Radians).
    double c = 2.0 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1.0 - a));

    // Distance.
    // const Double kEarthRadiusMiles = 3956.0;
    const Double kEarthRadiusKms = 6376.5;
    dDistance = kEarthRadiusKms * c;

    return dDistance;
}

如果距离矩形任意顶点之间的距离小于圆的半径则圆和矩形重叠。 如果圆心和所有顶点之间的距离大于圆的半径,并且所有这些距离都小于矩形的宽度和高度,则圆应该在矩形的内部。

如果您发现问题,请随时更正我的代码,因为我确信存在一些我没有想到的情况。

另外,我不确定这是否适用于跨越半球两端的矩形,因为距离方程可能会崩溃。

public string Test(double cLat,
    double cLon,
    double cRadius,
    double rlat1,
    double rlon1,
    double rlat2,
    double rlon2,
    double rlat3,
    double rlon3,
    double rlat4,
    double rlon4)
{
    double d1 = Calc(cLat, cLon, rlat1, rlon1);
    double d2 = Calc(cLat, cLon, rlat2, rlon2);
    double d3 = Calc(cLat, cLon, rlat3, rlon3);
    double d4 = Calc(cLat, cLon, rlat4, rlon4);

    if (d1 <= cRadius ||
        d2 <= cRadius ||
        d3 <= cRadius ||
        d4 <= cRadius)
    {

        return "Circle and Rectangle intersect...";
    }

    double width = Calc(rlat1, rlon1, rlat2, rlon2);
    double height = Calc(rlat1, rlon1, rlat4, rlon4);

    if (d1 >= cRadius &&
        d2 >= cRadius &&
        d3 >= cRadius &&
        d4 >= cRadius &&
        width >= d1 &&
        width >= d2 &&
        width >= d3 &&
        width >= d4 &&
        height >= d1 &&
        height >= d2 &&
        height >= d3 &&
        height >= d4)
    {
        return "Circle is Inside of Rectangle!";
    }



    return "NO!";
}

This should work for any points on earth. If you want to change it to a different size sphere just change the kEarchRadiusKms to whatever radius you want for your sphere.

This method is used to calculate the distance between to lat and lon points.

I got this distance formula from here:
http://www.codeproject.com/csharp/distancebetweenlocations.asp

public static double Calc(double Lat1, double Long1, double Lat2, double Long2)
{
    double dDistance = Double.MinValue;
    double dLat1InRad = Lat1 * (Math.PI / 180.0);
    double dLong1InRad = Long1 * (Math.PI / 180.0);
    double dLat2InRad = Lat2 * (Math.PI / 180.0);
    double dLong2InRad = Long2 * (Math.PI / 180.0);

    double dLongitude = dLong2InRad - dLong1InRad;
    double dLatitude = dLat2InRad - dLat1InRad;

    // Intermediate result a.
    double a = Math.Pow(Math.Sin(dLatitude / 2.0), 2.0) +
               Math.Cos(dLat1InRad) * Math.Cos(dLat2InRad) *
               Math.Pow(Math.Sin(dLongitude / 2.0), 2.0);

    // Intermediate result c (great circle distance in Radians).
    double c = 2.0 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1.0 - a));

    // Distance.
    // const Double kEarthRadiusMiles = 3956.0;
    const Double kEarthRadiusKms = 6376.5;
    dDistance = kEarthRadiusKms * c;

    return dDistance;
}

If the distance between any vertex of the rectangle is less than the distance of the radius of the circle then the circle and rectangle overlap. If the distance between the center of the circle and all of the vertices is greater than the radius of the circle and all of those distances are shorter than the width and height of the rectangle then the circle should be inside of the rectangle.

Feel free to correct my code if you can find a problem with it as I'm sure there some condition that I have not thought of.

Also I'm not sure if this works for a rectangle that spans the ends of the hemispheres as the distance equation might break down.

public string Test(double cLat,
    double cLon,
    double cRadius,
    double rlat1,
    double rlon1,
    double rlat2,
    double rlon2,
    double rlat3,
    double rlon3,
    double rlat4,
    double rlon4)
{
    double d1 = Calc(cLat, cLon, rlat1, rlon1);
    double d2 = Calc(cLat, cLon, rlat2, rlon2);
    double d3 = Calc(cLat, cLon, rlat3, rlon3);
    double d4 = Calc(cLat, cLon, rlat4, rlon4);

    if (d1 <= cRadius ||
        d2 <= cRadius ||
        d3 <= cRadius ||
        d4 <= cRadius)
    {

        return "Circle and Rectangle intersect...";
    }

    double width = Calc(rlat1, rlon1, rlat2, rlon2);
    double height = Calc(rlat1, rlon1, rlat4, rlon4);

    if (d1 >= cRadius &&
        d2 >= cRadius &&
        d3 >= cRadius &&
        d4 >= cRadius &&
        width >= d1 &&
        width >= d2 &&
        width >= d3 &&
        width >= d4 &&
        height >= d1 &&
        height >= d2 &&
        height >= d3 &&
        height >= d4)
    {
        return "Circle is Inside of Rectangle!";
    }



    return "NO!";
}
最近可好 2024-07-17 16:07:41

再尝试一次......

我认为解决方案是测试一组点,正如 Jason S 所建议的那样,但我不同意他选择的点,我认为这在数学上是错误的。

您需要找到纬度/经度框两侧到圆心的距离为局部最小值或最大值的点。 将这些点添加到角点集中,然后上面的算法应该是正确的。

即,让经度为 x 维度,纬度为 y 维度,让每个维度
盒子的一侧是参数曲线 P(t) = P0 + t (P1-P0) for o <= t <= 1.0,其中
P0和P1是两个相邻的角。

令 f(P) = f(Px, Py) 为距圆心的距离。

那么f(P0+t(P1-P0))就是t的距离函数:g(t)。 找到距离函数的导数为零的所有点:g'(t) == 0。(当然,丢弃解决方案超出了域 0 <= t <= 1.0)不幸的是

,这需要找到零的先验表达式,因此没有封闭形式的解决方案。 此类方程只能通过牛顿-拉夫逊迭代来求解。

好吧,我知道你想要的是代码,而不是数学。 但数学是我所拥有的一切。

One more try at this...

I think the solution is to test a set of points, just as Jason S has suggested, but I disagree with his selection of points, which I think is mathematically wrong.

You need to find the points on the sides of the lat/long box where the distance to the center of the circle is a local minimum or maximum. Add those points to the set of corners and then the algorithm above should be correct.

I.e, letting longitude be the x dimension and latitude be the y dimension, let each
side of the box be a parametric curve P(t) = P0 + t (P1-P0) for o <= t <= 1.0, where
P0 and P1 are two adjacent corners.

Let f(P) = f(P.x, P.y) be the distance from the center of the circle.

Then f (P0 + t (P1-P0)) is a distance function of t: g(t). Find all the points where the derivative of the distance function is zero: g'(t) == 0. (Discarding solutions outsize the domain 0 <= t <= 1.0, of course)

Unfortunately, this needs to find the zero of a transcendental expression, so there's no closed form solution. This type of equation can only solved by Newton-Raphson iteration.

OK, I realize that you wanted code, not the math. But the math is all I've got.

腻橙味 2024-07-17 16:07:41

有关欧几里得几何答案,请参阅:圆-矩形碰撞检测(交集)

For the Euclidean geometry answer, see: Circle-Rectangle collision detection (intersection)

~没有更多了~
我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
原文