How to find Bounding Box given 4 Lat/Long and Bearing? [duplicate]

Question

I need to find a bounding box given 4 lat/long points and a bearing (as seen in the example picture). I always know which two points are lined up by the bearing (1 and 2 in the example), so I will always know the length of the bounding box. The width however is arbitrary, with the points being anywhere along the lines (3 and 4 in the example).

My first thought is that I'll have to calculate the angles between the points (1 & 3, 1 & 4, 2 & 3, 2 & 4) and then use a series of "law of cosine" equations to calculate the corner points. Is there a simpler way? Would that even work?

Answer 1

If you can offer an azimuth, you could generate lines from the points and azimuth. then intersect to find corners. Once corners are generated, build the bounding box from these points to include points 1-4 Ex: point 1 azimuth generates line A point 2 azimuth generates line B point 3 azimuth generates line C point 4 azimuth generates line D

intersect lines A:C A:D D:B C:D Create points from the intersections

If you're using QGIS: Under research tools > Polygon From Extent InputLayer: Your points select a location for an output polygon shapefile...

the output will be a bounding box that extends no further then the points features' extent... while containing all the point features.

UPDATE: Here is that screenshot I promised, that is creating a Polygon from Extent using QGIS

Answer 2

The question makes little or no sense except in projected coordinates, because a "bounding box" and even a constant bearing have ambiguous meanings on a spheroid. In the projected plane, all you have to do is rotate the points to make the bearing horizontal or vertical. Finding the bounding box is then done the usual way (by obtaining the extreme values of the coordinates). Rotating this figure back finishes the job.

To maintain high precision, rotate the points around some central (or nearby) location.

---

If we were to imagine a GIS that supported the basic operations of (a) rotating figures and (b) finding bounding boxes--which most of them do--its solution would read like this:

BoundingBox = Rotate(Extent(Rotate(points, bearing)), -bearing)

where

Rotate(p, a)

rotates a feature p by an amount a and

Extent(p)

returns the (rectangular) extent of a feature.

---

Here is an example showing the rotated and original situations:

R code (which is readily ported to Python or any other platform supporting basic matrix operations) follows. Most of it just generates sample data and plots the results.

#
# Sample data (looking like those of the question,
# with typical ranges of projected coordinates).
#
xy <- cbind(1:4 + 200000, c(0,5,0,4) + 4000000)
bearing <- 45 * pi/180 # Radians east of north
#
# Helper functions.
#
rotate <- function(xy, a, origin=c(0,0)) {
  c <- cos(a); s <- sin(a)
  return(t(matrix(c(c, -s, s, c), 2) %*% (t(xy) - origin) + origin))
}
extent <- function(xy) {
  e <- apply(xy, 2, range)
  return(matrix(c(e[1,1], e[2,1], e[2,1], e[1,1],
                  e[1,2], e[1,2], e[2,2], e[2,2]), ncol=2, byrow=FALSE))
}
#
# Compute the oriented bounding box.
#
center <- apply(xy, 2, mean)
bb <- rotate(extent(rotate(xy, bearing, center)), -bearing, center)
#
# Display the points and their oriented bounding box.
#
plot(rbind(bb, xy), type="n", asp=1, xlab="X", ylab="Y", main="Solution")
polygon(bb, col="#f0f0f0")
points(xy, pch=19, col="Red")

Answer 3

@LMHall's answer pointed me in the right direction and I found a solution based on something by Chris Veness (here). Basically Chris published some scripts related to latitude/longitude calculations, one of which was how to calculate the intersection point given two points and their bearings. So to get the corners of the bounding box I just take each combination of top/bottom and left/right (1 & 3, 1 & 4, 2 & 3, 2 & 4) and find the intersections using the known bearing and adjusting accordingly. For example to find the bottom right of the image I'd calculate the intersection of points 1 & 3 using the bearing + 90 for the direction of point 1 and the bearing - 180 for the direction of point 3.

I can take no credit for the algorithm or even really explain it in terms of how it works geometrically, but its worked in my testing. Below is my java translation from the javascript version provided by Chris

public static CoordD getIntersection(CoordD point1, double bearing1, CoordD point2, double bearning2) {
    double lat1 = rad(point1.latitude); double lon1 = rad(point1.longitude);
    double lat2 = rad(point2.latitude); double lon2 = rad(point2.longitude);
    double bearing13 = rad(bearing1); double bearing 23 = rad(bearing2);
    double dLat = lat2 - lat1; double dLon = lon2 - lon1;

    double dist12 = 2 * Math.asin( Math.sqrt( Math.sin(dLat / 2) * Math.sin(dLat / 2) +
        Math.cos(lat1) * Math.cos(lat2) * Math.sin(dLon / 2) * Math.sin(dLon / 2) ) );
    if (dist12 == 0) return null;

    double bearingA = Math.acos( ( Math.sin(lat2) - Math.sin(lat1) * Math.cos(dist12) ) /
        ( Math.sin(dist12) * Math.cos(lat1) ) );
    double bearingB = Math.acos( ( Math.sin(lat1) - Math.sin(lat2) * Math.cos(dist12) ) /
        ( Math.sin(dist12) * Math.cos(lat2) ) );
    if (Double.isNaN(bearingA)) bearingA = 0;
    if (Double.isNaN(bearingB)) bearingB = 0;

    double bearing12, bearing21;
    if (Math.sin(dLon) > 0) {
        bearing12 = bearingA;
        bearing21 = 2 * Math.PI - bearingB;
    } else { 
        bearing12 = 2 * Math.PI - bearingA;
        bearing21 = bearingB;
    }

    double alpha1 = (bearing13 - bearing12 + Math.PI) % (2 * Math.PI) - Math.PI; // Angle 2-1-3
    double alpha2 = (bearing21 - bearing23 + Math.PI) % (2 * Math.PI) - Math.PI; // Angle 1-2-3

    if (Math.sin(alpha1) == 0 && Math.sin(alpha2) == 0) return null; // Infinite intersections
    if (Math.sin(alpha1) * Math.sin(alpha2) < 0) return null; // Ambiguous intersection

    // needed?
    // alpha1 = Math.abs(alpha1);
    // alpha2 = Math.abs(alpha2);

    double alpha3 = Math.acos( -Math.cos(alpha1) * Math.cos(alpha2) +
        Math.sin(alpha1) * Math.sin(alpha2) * Math.cos(dist12) );
    double dist13 = Math.atan2( Math.sin(dist12) * Math.sin(alpha1) * Math.sin(alpha2),
        Math.cos(alpha2) + Math.cos(alpha1) * Math.cos(alpha3) );

    double lat3 = Math.asin( Math.sin(lat1) * Math.cos(dist13) +
        Math.cos(lat1) * Math.sin(dist13) * Math.cos(bearing13) );

    double dLon13 = Math.atan2( Math.sin(bearing13) * Math.sin(dist13) * Math.cos(lat1),
        Math.cos(dist13) - Math.sin(lat1) * Math.sin(lat3) );
    double lon3 = lon1 + dLon3;
    lon3 = (lon3 + 3 * Math.PI) % ( 2* Math.PI) - Math.PI // normalize to +/-180

    return new CoordD(deg(lat3), deg(lon3));
}

rad() and deg() are just helper functions that translate between radians and degrees. CoordD is a helper class that just contains two double to store a lat/long point.