Georreferenciación directa con X,Y,Z y Omega, Phi, Kappa

Question

Tengo imágenes tomadas desde el avión, es longitud, latitud, elevación, omega, phi, kappa. Supongo que la tierra es plana (no necesito ningún DEM en ese escenario) y quiero calcular cómo se clava la imagen en el suelo. Quiero encontrar el centro y las 4 esquinas de la imagen clavada en el suelo.

¿Podría aconsejar alguna fórmula sencilla de cálculo? O una herramienta gratuita que tenga esos cálculos para comprobar los resultados de mi propia implementación. He intentado utilizar GRASS pero no he tenido suerte, es demasiado complicado :(

Gracias de antemano.

ACTUALIZACIÓN: Las imágenes se tomaron desde unos 3,2k metros sin GCPs físicos. La precisión que obtuvimos es de unos 20 metros con ángulos de cámara (Omega, Phi) de unos 2 grados. ¿Cree usted que es posible obtener una mejor precisión? Estamos usando C# y aquí está nuestra implementación:

    /// <summary>
    /// 
    /// </summary>
    /// <param name="omega">camera angle, in grads</param>
    /// <param name="phi">camera angle, in grads</param>
    /// <param name="kappa">camera angle, in grads</param>
    /// <param name="longitude">camera position</param>
    /// <param name="latitude">camera position</param>
    /// <param name="altitude">camera position</param>
    /// <param name="altOverGround">average flight height, it is defined on flight planning stage</param>
    /// <param name="imageLength">it is defined on flight planning stage, calculated from altOverGround and focal lenght, it is in meters</param>
    /// <param name="imageWidth">it is defined on flight planning stage, calculated from altOverGround and focal lenght, it is in meters</param>
    /// <returns></returns>
    public static DbGeography CalculateProjection(double omega, double phi, double kappa, double longitude, double latitude, double altitude,
        double altOverGround, double imageLength, double imageWidth)
    {
        // coordinates of the camera
        DbGeography cameraCoordinates = CreatePoint(longitude, latitude);

        // converting angles from grad to radian
        omega = ConvertFromGradToRadian(omega);
        phi = ConvertFromGradToRadian(phi);
        kappa = ConvertFromGradToRadian(kappa);

        // get image center (projected to groud) as offset to camera coordinates
        // result of imageCenter is in image coordinate system
        Point imageCenter = GetImageCenter(omega, phi, kappa, altitude);

        // calculating real image length for this height
        var realImageLength = imageCenter.Z * imageLength / altOverGround;
        var realImageWidth = imageCenter.Z * imageWidth / altOverGround;

        // corners offset in image coordinates system
        var xOfsfet = realImageLength / 2;
        var yOffset = realImageWidth / 2;

        // corners list
        List<Point> imageCorners = new List<Point>() { 
            new Point(xOfsfet, yOffset, imageCenter.Z),
            new Point(-xOfsfet, yOffset, imageCenter.Z),
            new Point(-xOfsfet, -yOffset, imageCenter.Z),
            new Point(xOfsfet, -yOffset, imageCenter.Z)
        };

        List<DbGeography> projectionCorners = new List<DbGeography>();
        foreach (var point in imageCorners)
        {
            // get image corner projected to the ground
            Point projectedPoint = GetProjectionCoordinates(omega, phi, kappa, altitude, point);

            // add projected point to list of projected corners
            projectionCorners.Add(CreatePointWithOffset(cameraCoordinates, projectedPoint.Y, projectedPoint.X));
        }

        // add first point to complete the polygon
        projectionCorners.Add(projectionCorners[0]);

        return CreatePolygonFromListOfPoints(projectionCorners);            
    }

    private static Point GetProjectionCoordinates(double omega, double phi, double kappa, double altitude, Point imageCorner)
    {
        // converting point coordinates from image's coordinate system to normal coordinate system
        Point imagePoint = ConvertImageCoordToNormal(omega, phi, kappa, imageCorner);

        // get coordinates of the projected point
        return GetProjectionCoordinates(altitude, imagePoint);
    }

    private static Point GetProjectionCoordinates(double altitude, Point imagePoint)
    {
        var x = altitude * imagePoint.X / imagePoint.Z;
        var y = altitude * imagePoint.Y / imagePoint.Z;
        return new Point(x, y, altitude);
    }

    private static Point ConvertImageCoordToNormal(double omega, double phi, double kappa, Point imageCorner)
    {
        double x1, x2, y1, y2, z1, z2, x, y, z;
        kappa = Math.PI / 2 - kappa;
        x1 = imageCorner.X * Math.Cos(kappa) + imageCorner.Y * Math.Sin(kappa);
        y1 = -imageCorner.X * Math.Sin(kappa) + imageCorner.Y * Math.Cos(kappa);
        z1 = imageCorner.Z;

        x2 = x1 * Math.Cos(phi) - z1 * Math.Sin(phi);
        y2 = y1;
        z2 = x1 * Math.Sin(phi) + z1 * Math.Cos(phi);

        x = x2;
        y = y2 * Math.Cos(omega) + z2 * Math.Sin(omega);
        z = -y2 * Math.Sin(omega) + z2 * Math.Cos(omega);
        return new Point(x, y, z);
    }

    /// <summary>
    /// Calculates altitude of the image center
    /// </summary>
    /// <param name="omega"></param>
    /// <param name="phi"></param>
    /// <param name="kappa"></param>
    /// <param name="altitude"></param>
    /// <returns>Altitude of image center</returns>
    private static Point GetImageCenter(double omega, double phi, double kappa, double altitude)
    {
        // calculating coordinates of image center
        var xCenter = altitude * Math.Tan(phi);
        var yCenter = altitude * Math.Tan(omega);
        var imageAltitude = altitude / (Math.Cos(omega) * Math.Cos(phi));

        // recalculating center with kappa value
        var alpha = Math.Atan(yCenter / xCenter) + kappa;
        var aa = Math.Sqrt(xCenter * xCenter + yCenter * yCenter);

        double imageCenterX = Math.Cos(alpha) * aa;
        double imageCenterY = Math.Sin(alpha) * aa;
        return new Point(imageCenterX, imageCenterY, imageAltitude);
    }

    private static double ConvertFromGradToRadian(double omega)
    {
        omega = Math.PI * omega / 200;
        return omega;
    }

y la clase de punto:

class Point
{
    public double X { get; set; }
    public double Y { get; set; }
    public double Z { get; set; }
    public Point(double x, double y, double z)
    {
        X = x;
        Y = y;
        Z = z;
    }
}

Answer 1

Aquí está La forma más fácil y gratuita de georreferenciar un archivo de imagen utilizando QGIS y ici son más enlaces

Georreferenciar una imagen en QGis es una tarea muy sencilla. A continuación se explica cómo hacerlo:

Ejecute QGis y cree un nuevo proyecto en blanco. No necesita guardarlo. Deberías tener ya algún otro material georreferenciado, lo mejor son las fotos aéreas, pero también las formas podrían servir. Abre la carpeta donde están esos archivos, y suéltalos en el proyecto en blanco. QGis cargará y los mostrará. En la ventana de QGis, haga clic en Plugins - Georeferenciador - Georeferenciador. La herramienta de georreferenciación se iniciará y mostrará una ventana. Haga clic en Archivo - Abrir ráster, y busque el archivo archivo de imagen, en mi caso un archivo jpg. El plugin carga el archivo en su ventana. Haga clic en la herramienta Añadir punto de la barra de herramientas (es el icono con tres puntos rojos), y luego haga clic en un punto conocido en el mapa. Es debe ser un punto que pueda reconocer en su otro mapa, como el borde de una casa, o una característica natural especial, etc. Aparece una ventana emergente que le pedirá que introduzca las coordenadas. Haga clic en el botón "Desde el mapa". La ventana emergente se cierra temporalmente y se vuelve al mapa original. original. Navegue con el teclado y la rueda del ratón y localice el punto en el mapa correspondiente al punto sobre el que hizo clic en el punto anterior punto. Haga clic allí y volverá a la ventana emergente. Un punto sobre el punto que acaba de georreferenciar. Repita los los puntos 5 y 6 en varios puntos diferentes. Debe elegir puntos en diferentes zonas del mapa, y no muy próximos entre sí. Cuanto más distancia, mejor. Son necesarios al menos tres. Una vez que hayas terminado haga clic en Archivo - Iniciar georreferenciación. Aparece una ventana emergente. Deberá rellenar al menos el archivo "Output raster". Se generará un archivo tiff georeferenciado. También puede jugar con la casilla "Transformación tipo". A mí, por ejemplo, me funcionó mejor la transformación Helmert. Haz clic en OK y el plugin procesará tu archivo. Ahora deberías cerrar la ventana (puede guardar sus puntos en un archivo separado, en en un archivo separado, en caso de que usted necesita para repetir el proceso más tarde, que recomiendo hacer). Por último, debes comprobar tu mapa. Una forma rápida y sencilla es utilizar la ventana principal, que ya contiene su mapa de referencia. Ahora, arrastre el archivo tiff que generó en la ventana principal: el mapa se añadido. Haz doble clic en su disposición en el panel izquierdo, y mueve el botón Transparencia global al 50%. Ahora puede comparar el mapa original y el nuevo, para ver si coinciden. Si Si no es así, puede volver a los puntos anteriores y ajustar algunos parámetros. Su archivo tiff está ahora georreferenciado, utilizando el mismo sistema de coordenadas que el mapa que utilizó como referencia. Ahora está listo para ser importado en mapserver.