Home » Software » PostGIS: Spatially enabled Relational Database Sytem » Affine Transformation Operations in PostGIS » Transformation parameters query

Transformation parameters query

The following is a proof of concept query showing how a PostgreSQL query could give the transformation parameters required for an affine transformation. The source of the procedure is from GRASS at transform.c.

Note that the v.transform code expects the input file (control points) to be in the form ax ay bx by, where 'a' is the starting coordinate system (the bad coordinates in the previous example) and 'b' is the target system (the good coordinates). The

SQL version is as follows:

Sample Session

1.

CREATE TABLE link (
    gid integer NOT NULL,
    a_x double precision,
    a_y double precision,
    b_x double precision,
    b_y double precision
);
-- load data into table
-- see attached file at bottom of page
-- possibly use the file 'grass_control_pts.csv' ???
COPY link FROM 'control_points.txt' ;

2.

CREATE TEMPORARY TABLE cc AS
SELECT
(SELECT count(b_x) FROM link) AS cc00,
(SELECT sum(b_x) FROM link) AS cc01,
(SELECT sum(b_y) FROM link) AS cc02,
sum(l1.b_x * l2.b_x) AS cc11,
sum(l1.b_x * l2.b_y) AS cc12,
sum(l1.b_y * l2.b_y) AS cc22
 FROM link l1 , link l2 WHERE l1.gid = l2.gid;

3.

CREATE TEMPORARY TABLE aa AS
SELECT
(SELECT sum(a_y) FROM link) AS aa0,
sum(l1.a_y * l2.b_x) AS aa1,
sum(l1.a_y * l2.b_y) AS aa2
 FROM link l1 , link l2 WHERE l1.gid = l2.gid;

4.

CREATE TEMPORARY TABLE bb AS
SELECT
(SELECT sum(a_x) FROM link) AS bb0,
sum(l1.a_x * l2.b_x) AS bb1,
sum(l1.a_x * l2.b_y) AS bb2
 FROM link l1 , link l2 WHERE l1.gid = l2.gid;

5.

CREATE TEMPORARY TABLE det AS
SELECT
cc00*cc11*cc22 +
cc01*cc12*cc02 +
cc02*cc01*cc12 -
cc00*cc12*cc12 -
cc01*cc01*cc22 -
cc02*cc11*cc02 AS det FROM cc;

6.

CREATE TEMPORARY TABLE inv_cc AS
SELECT
(cc11*cc22-cc12*cc12)/det AS cc00,
(cc12*cc02-cc01*cc22)/det AS cc01,
(cc01*cc12-cc11*cc02)/det AS cc02,
(cc00*cc22-cc02*cc02)/det AS cc11,
(cc01*cc02-cc00*cc12)/det AS cc12,
(cc00*cc11-cc01*cc01)/det AS cc22
FROM cc,det;

7.

SELECT
(aa0*cc00+aa1*cc01+aa2*cc02) AS b0,
(aa0*cc01+aa1*cc11+aa2*cc12) AS b1,
(aa0*cc02+aa1*cc12+aa2*cc22) AS b2,

(bb0*cc00+bb1*cc01+bb2*cc02) AS b3,
(bb0*cc01+bb1*cc11+bb2*cc12) AS b4,
(bb0*cc02+bb1*cc12+bb2*cc22) AS b5
FROM aa,bb,inv_cc;

This query requires a table called link with the following fields - gid (primary key), a_x, a_y , b_x, b_y. The 'a' values are the 'from' coordinates and the 'b' values are the 'to' coordinates. Using the attached control points the result of this query is:

        b0        |         b1          |        b2         |        b3        |        b4        |         b5          
------------------+---------------------+-------------------+------------------+------------------+---------------------
 -28138.394850347 | -0.0135202854235867 | 0.997400773420259 | 5017.08164289594 | 1.00231638907948 | 0.00918961946271679

These are the exact results from 'R'

Attachments:

control_points.txt
grass_control_pts.csv_.txt