Calculating coefficients for Affine Transformation

Introduction

Geospatial rasters inherently have two coordinate systems associated with them: pixel indices (i,j) and real world coordinates (x,y). For the simple case covered on this page, where a linear transform is sufficient to relate the two coordinate systems, physically significant parameters may be used to define the linear relation. These physically significant parameters control the magnitude and orientation of the basis vectors i_b and j_b defined in the following figure:

In the above image, the x and y axes represent geocoordinates (e.g., easting and northing in a projected system or longitude and latitude in a geographic system.) The basis vectors i_b and j_b are the transformed unit vectors i and j, respectively. The magnitude of the basis vectors represent the distance between grid cells in the units of the geocoordinate system. The angle between the x axis and basis vector i_b is given by θ_i. The angle between the two basis vectors is designated by θ_ij. In a "normal" right handed coordinate system, θ_ij is +90 degrees, but the value -90 degrees (indicating that the j axis needs to be flipped) is also very common. Any value other than plus or minus 90 degrees indicates that the linear transform includes shearing parallel to the i axis.

These simple linear relationships between pixel coordinates and geocoordinates are typically represented as modular combinations of individual transformations. Regardless of the order in which they are combined, a linear transform results. The transform is then used to convert coordinates between the two coordinate systems of the raster. This page describes a set of ubiquitous individual transformations and combines them to produce a specific transformation which respects the physically significant parameters which a typical user is likely to be interested in.

The physically significant inputs which drive the model are:

• The magnitude of the i_b basis vector (distance between pixels along i axis)
• The magnitude of the j_b basis vector (distance between pixels along j axis)
• The angle θ_i (the angle by which the raster's grid needs to be rotated, positive clockwise -- for compatibility with "heading/bearing")
• The angle θ_ij (the angle from the i basis vector to the j basis vector positive counterclockwise -- for consistency with right-handed coordinate system)

The level of math required to follow along is advanced high school algebra or introductory college algebra. It should be accessible to anyone with a science, math, or engineering background. Lacking this, the nontechnical introduction to ​matrix multiplication on wikipedia should provide sufficient background.

Individual operations

This page discusses operations in two dimensions only. Each operation is presented as a 2x2 matrix, and each operation performs only one function. These operations were taken from ​wikipedia. While the matrices presented here contain the bulk of the functionality of a finished affine transform, they are not complete affine transforms themselves. Each transformation will be presented in matrix and equation form: they are equivalent representations.

Rotation

There are two directions one can rotate in two dimensions: clockwise and counter clockwise. The transformations are different. The counter-clockwise rotation is:

• x' = xcosθ − ysinθ
• y' = xsinθ + ycosθ

To rotate in the clockwise direction:

• x' = xcosθ + ysinθ
• y' = − xsinθ + ycosθ

Scaling

Scaling is used to set the size of the raster's grid cells in the x and y direction. The transformation is as follows:

• x' = x s_x
• y' = y s_y

Shearing

Shearing is visually equivalent to a "slanting" which is parallel to either the x or the y axis. This is a less common operation than rotation and scaling. These are presented as individual operations: one for each axis.

Shearing parallel to the x axis takes the following form.

• x' = x + k_xy
• y' = y

While shearing parallel to the y axis has this form:

• x' = x
• y' = k_yx + y

Reflection

Reflection across the x axis (or "flipping" the y axis) is accomplished with the following transform:

• x' = x
• y' = -y

Combining individual operations

Whenever more than one of the above operations is required, they may be combined using ​matrix multiplication. As an example, all of the above matrices will be combined into one. The result of such a combination is still not an affine transform, however. It is just a 2x2 matrix which has all the individual transformation functions aggregated into it.

We will be calculating a new matrix, O, which is the aggregate of the following individual operations:

1. Reflection across the i axis
2. Scaling along the i and j axes
3. Shearing parallel to the i axis
4. clockwise rotation around the origin (by θ_i)

We do this by multiplying the 2x2 matrices of the individual operations together, as follows:

We will perform the matrix multiplications on the right hand side one at a time.

The above matrix equation is shorthand for four equations: one equation each for o₁₁, o₁₂, o₂₁ and o₂₂.

• o₁₁ = s_i cos(θ_i)
• o₁₂ = k_i s_j f cos(θ_i) + s_j f sin(θ_i)
• o₂₁ = -s_i sin(θ_i)
• o₂₂ = -k_i s_j f sin(θ_i) + s_j f cos(θ_i)

Notice that none of the coefficients in the O matrix may be said to represent pure scaling, rotation or shearing. Rather, they all have components of each of these operations factored in. The final matrix equation has some terms which need to be calculated. This will be performed in the next section using the real input parameters. Once these terms have been calculated, they may be plugged into the above equations to arrive at the coefficients which must be stored in the geotransform matrix.

Also notice that it is not necessary to compute this matrix every time one wants to convert between pixel indices and geographic location. The coefficients are computed once for the entire raster, and may be reused for every pixel calculation. You would use this aggregate matrix O exactly as you would use any of the individual matrices:

• x' = o₁₁ i + o₁₂ j
• y' = o₂₁ i + o₂₂ j

Constructing the affine transformation

The 2x2 matrix O is the upper-left hand corner of the affine transformation, which is the 3x3 matrix, A. The right hand column contains the offsets, or translation, of the raster in the x and y directions. The bottom row is always filled with the numbers 0, 0, 1. The result of this is that there are six parameters to an affine transformation which can actually change, as shown:

• a₁₁ = o₁₁ = s_x ( (1 + k_x k_y) cosθ + k_y sinθ )
• a₁₂ = o₁₂ = s_x ( k_x cosθ + sinθ )
• a₁₃ = t_x
• a₂₁ = o₂₁ = s_y ( -(1 + k_x k_y) sinθ + k_y cosθ )
• a₂₂ = o₂₂ = s_y ( - k_x sinθ + cosθ )
• a₂₃ = t_y

where:

• s_x : scale factor in x direction
• s_y : scale factor in y direction
• t_x : offset in x direction
• t_y : offset in y direction
• θ : angle of rotation clockwise around origin
• k_x : shearing parallel to x axis
• k_y : shearing parallel to y axis

It is these six parameters (a₁₁...a₂₃)which are typically stored within a geospatial image file format to record the conversion from pixel index to geolocation. The actual conversion is described as follows:

• E = a₁₁ i + a₁₂ j + a₁₃
• N = a₂₁ i + a₂₂ j + a₂₃

where E is easting, N is northing, i is pixel column and j is pixel row. The last row of the matrix equation is always ignored, as it boils down to 1=1. It is there to make a square matrix used to calculate the inverse operation.

Link to postgis raster

The six parameters of the affine transform are given the following names in postgis raster:

• ScaleX = a₁₁
• SkewX = a₁₂
• OffsetX = a₁₃ = t_x
• SkewY = a₂₁
• ScaleY = a₂₂
• OffsetY = a₂₃ = t_y

With the exception of OffsetX and OffsetY, the names are somewhat arbitrary for the general case.