[1]:

import sympy

Notes on `PCi_j` Formalism for Overlappograms#

From the PCi_j formalism, the conversion between intermediate world coordinates and pixel coordinates is given by (from Greisen and Calabretta, 2002),

\[x_i = s_iq_i\]

\[q_i = m_{ij}(p_j - r_j)\]

or in more explicit matrix notation,

\[X = sQ\]

\[Q = MP\]

In other words, we transform our pixel coordinates \(P\) to our intermediate world coordinates \(Q\) through the transformation \(M\). This matrix \(M\) is typically referred to as the PCi_j matrix, where \(i\) indexes the world coordinates and \(j\) the pixel coordinates. \(s\) is a scalar that represents the CDELT conversion factor between pixel units and world units.

The resulting PC_ij matrix for our overlappogram assuming that the dispersion and latitude world axes are aligned with the y-like pixel axis, :math:`p_2`,

\[\begin{split}D(\mu) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & -\mu \\ 0 & 0 & 1 \end{bmatrix}\end{split}\]

where \(\mu\) is the order of the dispersion.

In general, the dispersion axis may be oriented at some angle \(\gamma\) relative to \(p_2\). Thus, before applying \(D\), we rotate \(P\) by \(\gamma\) into \(P^\prime\) such that \(P^\prime=R(\gamma)P\) where \(R\) is the usual rotation matrix,

\[\begin{split}R(\theta) = \begin{bmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{bmatrix}\end{split}\]

Lastly, we apply a rotation to align \(P^\prime\) with the intermediate world axes \(Q\) such that \(Q=R(\beta)P^\prime\). Note that this angle \(\beta\) is the relative orientation between \(P^\prime\) and \(Q\). However, we will typically be measuring the satellite roll angle relative to the original pixel axes \(P\). Thus, \(\beta=\alpha-\gamma\) where \(\alpha\) is the orientation of solar north relative to \(p_1\).

The full PC matrix can thus be written as,

\[M = R(\alpha-\gamma)D(\mu)R(\gamma)\]

which can be written explicitly as,

\[\begin{split}M = \begin{bmatrix} \cos{\left(\alpha \right)} & -\sin{\left(\alpha \right)} & - \mu\cos{\left(\alpha - \gamma \right)}\\ \sin{\left(\alpha \right)} & \cos{\left(\alpha \right)} & - \mu\sin{\left(\alpha - \gamma \right)}\\ 0 & 0 & 1 \end{bmatrix}\end{split}\]

We can evaluate this matrix product with sympy and look at a few example cases.

[2]:

mu = sympy.symbols('mu')
gamma = sympy.symbols('gamma')
alpha = sympy.symbols('alpha')
theta = sympy.symbols('theta')
phi = sympy.symbols('phi')
rotation = sympy.Matrix([
    [sympy.cos(theta), -sympy.sin(theta), 0],
    [sympy.sin(theta), sympy.cos(theta), 0],
    [0, 0, 1]
])
dispersion = sympy.Matrix([
    [1, 0, -mu],
    [0, 1, 0],
    [0 ,0, 1],
])

[18]:

PC_matrix = rotation.subs(theta, (alpha-gamma)) * dispersion * rotation.subs(theta, gamma)
PC_matrix.simplify()
sympy.printing.pprint(PC_matrix)

⎡cos(α)  -sin(α)  -μ⋅cos(α - γ)⎤
⎢                              ⎥
⎢sin(α)  cos(α)   -μ⋅sin(α - γ)⎥
⎢                              ⎥
⎣  0        0           1      ⎦

We can evaluate the PCi_j matrix at a few different combinations of \(\alpha,\gamma,\mu\).

[4]:

for a,g,m in [(0,0,0),
              (0,0,1),
              (sympy.pi/2,0,1),
              (0,sympy.pi/2,1),
              (sympy.pi/2,sympy.pi/2,1),
              (sympy.pi/4,0,1),
              (0,sympy.pi/4,1),
              (sympy.pi/4,sympy.pi/4,1)]:
    print(f'alpha={a}, gamma={g}, mu={m}')
    print('-----------------------------')
    sympy.printing.pprint(PC_matrix.subs([(alpha, a),(gamma, g),(mu, m)]))

alpha=0, gamma=0, mu=0
-----------------------------
⎡1  0  0⎤
⎢       ⎥
⎢0  1  0⎥
⎢       ⎥
⎣0  0  1⎦
alpha=0, gamma=0, mu=1
-----------------------------
⎡1  0  -1⎤
⎢        ⎥
⎢0  1  0 ⎥
⎢        ⎥
⎣0  0  1 ⎦
alpha=pi/2, gamma=0, mu=1
-----------------------------
⎡0  -1  0 ⎤
⎢         ⎥
⎢1  0   -1⎥
⎢         ⎥
⎣0  0   1 ⎦
alpha=0, gamma=pi/2, mu=1
-----------------------------
⎡1  0  0⎤
⎢       ⎥
⎢0  1  1⎥
⎢       ⎥
⎣0  0  1⎦
alpha=pi/2, gamma=pi/2, mu=1
-----------------------------
⎡0  -1  -1⎤
⎢         ⎥
⎢1  0   0 ⎥
⎢         ⎥
⎣0  0   1 ⎦
alpha=pi/4, gamma=0, mu=1
-----------------------------
⎡√2  -√2   -√2 ⎤
⎢──  ────  ────⎥
⎢2    2     2  ⎥
⎢              ⎥
⎢√2   √2   -√2 ⎥
⎢──   ──   ────⎥
⎢2    2     2  ⎥
⎢              ⎥
⎣0    0     1  ⎦
alpha=0, gamma=pi/4, mu=1
-----------------------------
⎡      -√2 ⎤
⎢1  0  ────⎥
⎢       2  ⎥
⎢          ⎥
⎢       √2 ⎥
⎢0  1   ── ⎥
⎢       2  ⎥
⎢          ⎥
⎣0  0   1  ⎦
alpha=pi/4, gamma=pi/4, mu=1
-----------------------------
⎡√2  -√2     ⎤
⎢──  ────  -1⎥
⎢2    2      ⎥
⎢            ⎥
⎢√2   √2     ⎥
⎢──   ──   0 ⎥
⎢2    2      ⎥
⎢            ⎥
⎣0    0    1 ⎦

We can do some simple computations that shows how this transforms

[19]:

p = sympy.symarray('p', 3)
r = sympy.symarray('r', 3)
pix_coord = sympy.Matrix([
    p[0] - r[0],
    p[1] - r[1],
    p[2] - r[2],
])

[20]:

sympy.printing.pprint(pix_coord)

⎡p₀ - r₀⎤
⎢       ⎥
⎢p₁ - r₁⎥
⎢       ⎥
⎣p₂ - r₂⎦

[21]:

world_coord = PC_matrix * pix_coord

[22]:

sympy.printing.pprint(world_coord)

⎡-μ⋅(p₂ - r₂)⋅cos(α - γ) + (p₀ - r₀)⋅cos(α) - (p₁ - r₁)⋅sin(α)⎤
⎢                                                             ⎥
⎢-μ⋅(p₂ - r₂)⋅sin(α - γ) + (p₀ - r₀)⋅sin(α) + (p₁ - r₁)⋅cos(α)⎥
⎢                                                             ⎥
⎣                           p₂ - r₂                           ⎦

[10]:

data_shape = (1073, 2000, 750)

[25]:

Q = world_coord.subs([
    (alpha, 0),
    (gamma, 0),
    (mu, 1),
    (r[0], (data_shape[2] + 1)/2),
    (r[1], (data_shape[1] + 1)/2),
    (r[2], 1),#(data_shape[0] + 1)/2),
])
sympy.printing.pprint(Q)

⎡p₀ - p₂ - 374.5⎤
⎢               ⎥
⎢  p₁ - 1000.5  ⎥
⎢               ⎥
⎣    p₂ - 1     ⎦

[36]:

sympy.printing.pprint(Q.subs([
    (p[0], (data_shape[2] + 1)/2),
    (p[1], (data_shape[1] + 1)/2),
    (p[2], 100),
]))

⎡-99.0⎤
⎢     ⎥
⎢  0  ⎥
⎢     ⎥
⎣ 99  ⎦

[ ]:

Notes on PCi_j Formalism for Overlappograms#

Notes on `PCi_j` Formalism for Overlappograms#