[GRASS-dev] transform additional points from XY?

[Glynn Clements]

Hamish wrote:

> If a georeferencing points file and target are set up it is possible to
> use g.transform to say how good the fit is.
> 
> I'd like to go one further: I'd like to pass g.transform a handful of x,y
> pixels from the the same neighborhood as the GCPs and get an approximate
> projected easting and northing back.
> 
> If possible at all, I guess we have to write something new for
> g.transform; any good ideas on how to do this?

CRS_compute_georef_equations() calculates the forward and reverse
transformation coefficients, while CRS_georef() applies the
transformation to a single point.

E12/N12 are the coefficients for the forward transformation, while
E21/N21 are the coefficients for the reverse transformation.

CRS_compute_georef_equations() computes both pairs of arrays, while
CRS_georef() expects one pair to be supplied, depending upon whether
you want the forward or reverse transformation.

Either the E12/N12/E2/N21 arrays need to be moved out of
compute_transformation() and into global variables, or the code to
transform the points needs to be added to compute_transformation().

The transformations are:

order=1:

	e = [E0 E1][1].[1]
	    [E2  0][e] [n]
	
	n = [N0 N1][1].[1]
	    [N2  0][e] [n]

order=2:

	e = [E0 E1 E3][1 ] [1 ]
	    [E2 E4  0][e ].[n ]
	    [E5  0  0][e²] [n²]
	
	n = [N0 N1 N3][1 ] [1 ]
	    [N2 N4  0][e ].[n ]
	    [N5  0  0][e²] [n²]

order=3:

	e = [E0 E1 E3 E6][1 ] [1 ]
	    [E2 E4 E7  0][e ].[n ]
	    [E5 E8  0  0][e²] [n²]
	    [E9  0  0  0][e³] [n³]
	
	n = [N0 N1 N3 N6][1 ] [1 ]
	    [N2 N4 N7  0][e ].[n ]
	    [N5 N8  0  0][e²] [n²]
	    [N9  0  0  0][e³] [n³]

["." = dot-product, (AE).N = N'EA.]

IOW, order=1 and order=2 are equivalent to order=3 with the higher
coefficients equal to zero.

-- 
Glynn Clements <glynn at gclements.plus.com>