¿Cómo encontrar el cuadro delimitador teniendo en cuenta 4 Lat/Long y Bearing?

Question

Necesito encontrar un cuadro delimitador dados 4 puntos de latitud y longitud y un rumbo (como se ve en la imagen de ejemplo). Siempre sé qué dos puntos están alineados por el rumbo (1 y 2 en el ejemplo), por lo que siempre sabré la longitud del cuadro delimitador. La anchura, sin embargo, es arbitraria, ya que los puntos se encuentran en cualquier punto de las líneas (3 y 4 en el ejemplo).

Mi primera idea es que tendré que calcular los ángulos entre los puntos (1 y 3, 1 y 4, 2 y 3, 2 y 4) y luego utilizar una serie de ecuaciones de la "ley del coseno" para calcular los puntos de esquina. ¿Existe una forma más sencilla? ¿Funcionaría?

Answer 1

Si puede ofrecer un acimut, podría generar líneas a partir de los puntos y el acimut. luego se intersectan para encontrar las esquinas. Una vez generadas las esquinas, construya el cuadro delimitador a partir de estos puntos para incluir los puntos 1-4 Ej: el acimut del punto 1 genera la línea A el azimut del punto 2 genera la línea B el punto 3 acimut genera la línea C el punto 4 acimut genera la línea D

intersección de líneas A:C A:D D:B C:D Crea puntos a partir de las intersecciones

Si está utilizando QGIS: En herramientas de investigación > Polígono desde la extensión InputLayer: Sus puntos seleccione una ubicación para un shapefile de polígono de salida...

el resultado será un cuadro delimitador que no se extiende más allá de la extensión de las características de los puntos... mientras que contiene todas las características de los puntos.

ACTUALIZACIÓN: Aquí está la captura de pantalla que prometí, que es la creación de un polígono a partir de una extensión utilizando QGIS

Answer 2

La pregunta tiene poco o ningún sentido si no es en coordenadas proyectadas, porque un "cuadro delimitador" e incluso un rumbo constante tienen significados ambiguos en un esferoide. En el plano proyectado, sólo hay que girar los puntos para que el rodamiento sea horizontal o vertical. La búsqueda de la caja delimitadora se realiza entonces de la forma habitual (obteniendo los valores extremos de las coordenadas). La rotación de esta figura hacia atrás termina el trabajo.

Para mantener una alta precisión, gire los puntos alrededor de algún lugar central (o cercano).

---

Si imagináramos un SIG que soportara las operaciones básicas de (a) rotar figuras y (b) encontrar cajas delimitadoras -lo que hace la mayoría de ellos-, su solución sería la siguiente:

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

donde

Rotate(p, a)

gira una característica p por una cantidad a y

Extent(p)

devuelve la extensión (rectangular) de una característica.

---

Aquí hay un ejemplo que muestra las situaciones giradas y originales:

R (que es fácilmente trasladable a Python o a cualquier otra plataforma que soporte operaciones matriciales básicas). La mayor parte se limita a generar datos de muestra y a trazar los resultados.

#
# 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

La respuesta de @LMHall me indicó la dirección correcta y encontré una solución basada en algo de Chris Veness ( ici ). Básicamente, Chris publicó algunos scripts relacionados con los cálculos de latitud/longitud, uno de los cuales era cómo calcular el punto de intersección dados dos puntos y sus rumbos. Así que para obtener las esquinas del cuadro delimitador sólo tomo cada combinación de arriba/abajo e izquierda/derecha (1 y 3, 1 y 4, 2 y 3, 2 y 4) y encuentro las intersecciones usando el rumbo conocido y ajustando en consecuencia. Por ejemplo, para encontrar la parte inferior derecha de la imagen, calcularía la intersección de los puntos 1 y 3 utilizando el rumbo + 90 para la dirección del punto 1 y el rumbo - 180 para la dirección del punto 3.

No puedo atribuirme el mérito del algoritmo ni explicarlo realmente en términos de cómo funciona geométricamente, pero ha funcionado en mis pruebas. A continuación es mi traducción java de la versión javascript proporcionada por 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() y deg() son sólo funciones de ayuda que traducen entre radianes y grados. CoordD es una clase de ayuda que sólo contiene dos dobles para almacenar un punto lat/long.