3

After reading this post How to georeference a vector layer with control points? I decided to give it a try to move/rotation/scale a polygon using Postgis.

I have two polygons, polygon_from and polygon_to

Below are their 4 vertices' longtitude-latitude. I intend to move polygon_from to polygon_to position.

 gid |      from_x      |      from_y      |       to_x       |       to_y
-----+------------------+------------------+------------------+------------------
   1 | 103.771789705736 | 1.29874218723039 | 103.771893494233 | 1.29880598084823
   2 | 103.7718553679   | 1.2987732903611  | 103.771964782989 | 1.29879195095702
   3 | 103.771887622999 | 1.29870762819626 | 103.771951679355 | 1.29871997730869
   4 | 103.771820808866 | 1.29867710104945 | 103.771879806787 | 1.29873515522057

I copy and pasted a function (below) from this website http://casoilresource.lawr.ucdavis.edu/drupal/book/export/html/532

-- Function: trans_param(IN sql text, OUT a double precision, OUT b double precision, OUT d double precision, OUT e double precision, OUT xoff double precision, OUT yoff double precision)

-- DROP FUNCTION trans_param(IN sql text, OUT a double precision, OUT b double precision, OUT d double precision, OUT e double precision, OUT xoff double precision, OUT yoff double precision);

CREATE OR REPLACE FUNCTION trans_param(IN sql text, OUT a double precision, OUT b double precision, OUT d double precision, OUT e double precision, OUT xoff double precision, OUT yoff double precision) AS $BODY$ DECLARE
        cc_row record;

        cc_det float;

        inv_cc00 float;
        inv_cc01 float;
        inv_cc02 float;
        inv_cc11 float;
        inv_cc12 float;
        inv_cc22 float;

BEGIN
        EXECUTE 'SELECT 
        count(a.to_x) as cc00,
        sum(a.to_x) as cc01,
        sum(a.to_y) as cc02,
        sum(a.to_x * b.to_x) as cc11,
        sum(a.to_x * b.to_y) as cc12,
        sum(a.to_y * b.to_y) as cc22,
        sum(a.from_y) as aa0,
        sum(a.from_y * b.to_x) as aa1,
        sum(a.from_y * b.to_y) as aa2,
        sum(a.from_x) as bb0,
        sum(a.from_x * b.to_x) as bb1,
        sum(a.from_x * b.to_y) as bb2
        from (' || sql || ') a,(' || sql|| ') b WHERE a.gid = b.gid' INTO cc_row;


        SELECT INTO cc_det 
                cc_row.cc00*cc_row.cc11*cc_row.cc22 + 
                cc_row.cc01*cc_row.cc12*cc_row.cc02 + 
                cc_row.cc02*cc_row.cc01*cc_row.cc12 - 
                cc_row.cc00*cc_row.cc12*cc_row.cc12 - 
                cc_row.cc01*cc_row.cc01*cc_row.cc22 - 
                cc_row.cc02*cc_row.cc11*cc_row.cc02;

        SELECT INTO inv_cc00 (cc_row.cc11*cc_row.cc22-cc_row.cc12*cc_row.cc12)/cc_det;
        SELECT INTO inv_cc01 (cc_row.cc12*cc_row.cc02-cc_row.cc01*cc_row.cc22)/cc_det;
        SELECT INTO inv_cc02 (cc_row.cc01*cc_row.cc12-cc_row.cc11*cc_row.cc02)/cc_det;
        SELECT INTO inv_cc11 (cc_row.cc00*cc_row.cc22-cc_row.cc02*cc_row.cc02)/cc_det;
        SELECT INTO inv_cc12 (cc_row.cc01*cc_row.cc02-cc_row.cc00*cc_row.cc12)/cc_det;
        SELECT INTO inv_cc22 (cc_row.cc00*cc_row.cc11-cc_row.cc01*cc_row.cc01)/cc_det;

        SELECT INTO a cc_row.bb0*inv_cc01+cc_row.bb1*inv_cc11+cc_row.bb2*inv_cc12;
        SELECT INTO b cc_row.bb0*inv_cc02+cc_row.bb1*inv_cc12+cc_row.bb2*inv_cc22;
        SELECT INTO d cc_row.aa0*inv_cc01+cc_row.aa1*inv_cc11+cc_row.aa2*inv_cc12;
        SELECT INTO e cc_row.aa0*inv_cc02+cc_row.aa1*inv_cc12+cc_row.aa2*inv_cc22;
        SELECT INTO xoff cc_row.bb0*inv_cc00+cc_row.bb1*inv_cc01+cc_row.bb2*inv_cc02;
        SELECT INTO yoff cc_row.aa0*inv_cc00+cc_row.aa1*inv_cc01+cc_row.aa2*inv_cc02;

END;

$BODY$   LANGUAGE 'plpgsql' VOLATILE; ALTER FUNCTION trans_param(IN sql text, OUT a double precision, OUT b double precision, OUT d double precision, OUT e double precision, OUT xoff double precision, OUT yoff double precision) OWNER TO postgres;

SELECT trans_param('select * from above_table');

After applying the PostGIS function, I obtained a,b,d,e,x_off,y_off parameters as below

(0.0824142911101262,47.5134245217341,0.00103114836995566,0.594639499427103,-9.93126928180573,-0.124261944539057)

However when I apply those parameters to the polygon_from, it goes to a very strange lat/long position and becomes an irregular shape (not what I expected).

Update polygon_table Set the_geom=(Select st_affine(the_geom,0.0824142911101262,47.5134245217341,0.00103114836995566,0.594639499427103,-9.93126928180573,-0.124261944539057) Where gid=polygon_from_gid;

I am wondering how to use the above trans_param() function and those ST_Affine parameters to shift/rotate/scale a polygon?

Xianlin
  • 487
  • 7
  • 15
  • I had a look at http://casoilresource.lawr.ucdavis.edu/drupal/book/export/html/532 (apparently by someone named Bruce) where it compares transformations from Grass, R, PostGIS and his own plpgsql trans_param() function. It seems that those transfors each use different naming or ordering of parameters. Thus, he, you or I may have got confused somewhere! Some parameters may need re-ordering (but I cannot say which/how). – Martin F Nov 23 '13 at 02:31

1 Answers1

2

I had the same problem: the result of the trans_params function defined as we can read here https://casoilresource.lawr.ucdavis.edu/software/postgis-spatially-enabled-relational-database-sytem/affine-transformation-operations-postgis/user-defined-function-plpgsql-compute-transformation-parameters/ was incorrect and very different from the same result for example I obtained with the ogr2org command (with GPC parameter). I would apply that function to more than 50 control points in UTM metric coordinates. I discovered finally that the problem was due to the PRECISION of the variables declared in that function: just changing the type of variables from "float" to "numeric" (mainly for the "cc_det" variable) solved the problem.

I also simplified a bit the code, always taking in account the v.transform code from GRASS source code https://grass.osgeo.org/programming6/transform_8c_source.html, thanks to OPEN SOURCE !

finally the correct code is:

-- DROP FUNCTION public.trans_param(text);

CREATE OR REPLACE FUNCTION public.trans_param(
    IN sql text,
    OUT a numeric,
    OUT b numeric,
    OUT d numeric,
    OUT e numeric,
    OUT xoff numeric,
    OUT yoff numeric)
  RETURNS record AS
$BODY$
DECLARE
        cc_row record;

        cc_det numeric;

        inv_cc00 numeric;
        inv_cc01 numeric;
        inv_cc02 numeric;
        inv_cc11 numeric;
        inv_cc12 numeric;
        inv_cc22 numeric;

BEGIN
        EXECUTE 'SELECT
        count(a.to_x) as cc00,
        sum(a.to_x)::numeric as cc01,
        sum(a.to_y)::numeric as cc02,
        sum(a.to_x * a.to_x)::numeric as cc11,
        sum(a.to_x * a.to_y)::numeric as cc12,
        sum(a.to_y * a.to_y)::numeric as cc22,
        sum(a.from_y)::numeric as aa0,
        sum(a.from_y * a.to_x)::numeric as aa1,
        sum(a.from_y * a.to_y)::numeric as aa2,
        sum(a.from_x)::numeric as bb0,
        sum(a.from_x * a.to_x)::numeric as bb1,
        sum(a.from_x * a.to_y)::numeric as bb2
        from (' || sql || ') a' INTO cc_row;


        SELECT INTO cc_det
                (cc_row.cc00*cc_row.cc11*cc_row.cc22 +
                cc_row.cc01*cc_row.cc12*cc_row.cc02 +
                cc_row.cc02*cc_row.cc01*cc_row.cc12 -
                cc_row.cc00*cc_row.cc12*cc_row.cc12 -
                cc_row.cc01*cc_row.cc01*cc_row.cc22 -
                cc_row.cc02*cc_row.cc11*cc_row.cc02)::numeric;

        SELECT INTO inv_cc00 ((cc_row.cc11*cc_row.cc22-cc_row.cc12*cc_row.cc12)/cc_det)::numeric;
        SELECT INTO inv_cc01 ((cc_row.cc12*cc_row.cc02-cc_row.cc01*cc_row.cc22)/cc_det)::numeric;
        SELECT INTO inv_cc02 ((cc_row.cc01*cc_row.cc12-cc_row.cc11*cc_row.cc02)/cc_det)::numeric;
        SELECT INTO inv_cc11 ((cc_row.cc00*cc_row.cc22-cc_row.cc02*cc_row.cc02)/cc_det)::numeric;
        SELECT INTO inv_cc12 ((cc_row.cc01*cc_row.cc02-cc_row.cc00*cc_row.cc12)/cc_det)::numeric;
        SELECT INTO inv_cc22 ((cc_row.cc00*cc_row.cc11-cc_row.cc01*cc_row.cc01)/cc_det)::numeric;

        SELECT INTO a (cc_row.bb0*inv_cc01+cc_row.bb1*inv_cc11+cc_row.bb2*inv_cc12)::numeric;
        SELECT INTO b (cc_row.bb0*inv_cc02+cc_row.bb1*inv_cc12+cc_row.bb2*inv_cc22)::numeric;
        SELECT INTO d (cc_row.aa0*inv_cc01+cc_row.aa1*inv_cc11+cc_row.aa2*inv_cc12)::numeric;
        SELECT INTO e (cc_row.aa0*inv_cc02+cc_row.aa1*inv_cc12+cc_row.aa2*inv_cc22)::numeric;
        SELECT INTO xoff (cc_row.bb0*inv_cc00+cc_row.bb1*inv_cc01+cc_row.bb2*inv_cc02)::numeric;
        SELECT INTO yoff (cc_row.aa0*inv_cc00+cc_row.aa1*inv_cc01+cc_row.aa2*inv_cc02)::numeric;

END;

$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
francesco
  • 21
  • 2