Geostationary projection

From David's Wiki
\( \newcommand{\P}[]{\unicode{xB6}} \newcommand{\AA}[]{\unicode{x212B}} \newcommand{\empty}[]{\emptyset} \newcommand{\O}[]{\emptyset} \newcommand{\Alpha}[]{Α} \newcommand{\Beta}[]{Β} \newcommand{\Epsilon}[]{Ε} \newcommand{\Iota}[]{Ι} \newcommand{\Kappa}[]{Κ} \newcommand{\Rho}[]{Ρ} \newcommand{\Tau}[]{Τ} \newcommand{\Zeta}[]{Ζ} \newcommand{\Mu}[]{\unicode{x039C}} \newcommand{\Chi}[]{Χ} \newcommand{\Eta}[]{\unicode{x0397}} \newcommand{\Nu}[]{\unicode{x039D}} \newcommand{\Omicron}[]{\unicode{x039F}} \DeclareMathOperator{\sgn}{sgn} \def\oiint{\mathop{\vcenter{\mathchoice{\huge\unicode{x222F}\,}{\unicode{x222F}}{\unicode{x222F}}{\unicode{x222F}}}\,}\nolimits} \def\oiiint{\mathop{\vcenter{\mathchoice{\huge\unicode{x2230}\,}{\unicode{x2230}}{\unicode{x2230}}{\unicode{x2230}}}\,}\nolimits} \)

The Geostationary projection is used to project 2D satellite images, photographs of the globe back into spherical coordinates or equirectangular images.
Usually it is assumed the satellite is always above the equator and is oriented correctly such that the north-pole is up.

Projection

Input:
Phi, Theta spherical coordinates
Theta spherical coordinates of the satellite position
Output:
Angle (in x,y) where positive y corresponds to the north pole and negative y corresponds to the south pole.

private static void LatLngToAngles(float lat, float lng, NetCDFData projection, out int targetX, out int targetY)
{
    // Given a latitude, longitude pair
    // Calculate the cartesian coordinates of the point
    // y is the height from the equator
    // x is the distance from the center of the earth to the camera
    // z is left/right of the vector to the satellite (on the plane of the equator)
    // i.e. the camera (satellite) is at (H, 0, 0)
    // and the north pole is at (0, projection.semiMinorAxis,0)

    float longitudeDelta = (lng * Mathf.PI / 180) - (projection.lonOrigin * Mathf.PI / 180);
    longitudeDelta = (longitudeDelta + 9 * Mathf.PI) % (2 * Mathf.PI);
    longitudeDelta -= Mathf.PI;
    float latitudeInRad = lat * Mathf.PI / 180;
    if (Mathf.Abs(longitudeDelta) > Mathf.PI / 2)
    {
        targetX = -1;
        targetY = -1;
        return;
    }
    float yPos = projection.semiMinorAxis * Mathf.Sin(latitudeInRad);
    float xPos = projection.semiMajorAxis * Mathf.Cos(longitudeDelta) * Mathf.Cos(latitudeInRad);
    float zPos = projection.semiMajorAxis * Mathf.Sin(longitudeDelta) * Mathf.Cos(latitudeInRad);

    // Get the tangent plane at this point.
    // If the camera is on a particular side of the tangent plane then
    // the camera cannot see this point.
    Vector3 position = new Vector3(xPos, yPos, zPos);
    Vector3 normalVector = position.normalized;
    float dotProduct = Vector3.Dot(normalVector, new Vector3(projection.height, 0, 0) - position);
    if (dotProduct < 0)
    {
        targetX = -1;
        targetY = -1;
        return;
    }

    //Debug.Log("Ypos is " + yPos);
    //Debug.Log("Xpos is " + xPos);

    float angleY = Mathf.Atan(yPos / (projection.height - xPos));
    float angleX = Mathf.Atan(zPos / (projection.height - xPos));

    //Debug.Log("AngleY is " + angleY);

    targetX = (int)((angleX - projection.xOffset) / projection.xScaleFactor);
    targetY = (int)((angleY - projection.yOffset) / projection.yScaleFactor);

    //Debug.Log("TargetY is: " + targetY);
}


Inverse Projection

See Reference


Resources