Geoinquiets

CREATE USER datumtest LOGIN PASSWORD 'datumtest' NOSUPERUSER INHERIT NOCREATEDB NOCREATEROLE;
CREATE DATABASE datumtest OWNER datumtest;
\c datumtest
CREATE EXTENSION postgis;

CREATE TABLE ed50_test_points (id serial, geom geometry(Point, 23031));
INSERT INTO ed50_test_points (geom) VALUES
    (ST_GeomFromText('POINT(300000.000 4500000.000)', 23031)),
    (ST_GeomFromText('POINT(315000.000 4740000.000)', 23031)),
    (ST_GeomFromText('POINT(520000.000 4680000.000)', 23031)),
    (ST_GeomFromText('POINT(420000.000 4600000.000)', 23031));

CREATE TABLE etrs89_test_points (id serial, geom geometry(Point, 25831));
INSERT INTO etrs89_test_points (geom) VALUES
    (ST_GeomFromText('POINT(300000.000 4500000.000)', 25831)),
    (ST_GeomFromText('POINT(315000.000 4740000.000)', 25831)),
    (ST_GeomFromText('POINT(520000.000 4680000.000)', 25831)),
    (ST_GeomFromText('POINT(420000.000 4600000.000)', 25831));

ST_Affine(geom, a, b, d, e, xoff, yoff)

ST_Affine(geom, 1.0000015503712145, 0.000007587539798467343, -0.000007587539798467343, 1.0000015503712145, -129.549, -208.185)

-- Function: public.icgc_23031_to_25831(character varying, character varying, character varying)

-- DROP FUNCTION public.icgc_23031_to_25831(character varying, character varying, character varying);

CREATE OR REPLACE FUNCTION public.icgc_23031_to_25831(
    schema_name character varying,
    table_name character varying,
    column_name character varying)
  RETURNS text AS
$BODY$
DECLARE
	real_schema name;
	src_srid integer;
	dst_srid integer;
BEGIN
	-- Set SRID values
	src_srid := 23031;
	dst_srid := 25831;
	
	-- Find, check or fix schema_name
	IF ( schema_name != '' ) THEN
		real_schema = schema_name;
	ELSE
		SELECT INTO real_schema current_schema()::text;
	END IF;

	-- Check original SRID
	IF (SELECT count(*) = 0 FROM geometry_columns WHERE f_table_schema = real_schema AND f_table_name = table_name AND f_geometry_column = column_name AND srid = src_srid) THEN
		RAISE EXCEPTION 'table original SRID not matching required SRID %', src_srid;
		RETURN false;
	END IF;

	-- Set new SRID to table metadata
	PERFORM updategeometrysrid(schema_name, table_name, column_name, dst_srid);

	-- Reproject using affine transform
	EXECUTE 'UPDATE ' || quote_ident(real_schema) || '.' || quote_ident(table_name) || ' SET ' || quote_ident(column_name) || ' = ST_Affine('|| quote_ident(column_name) ||', 1.0000015503712145, 0.000007587539798467343, -0.000007587539798467343, 1.0000015503712145, -129.549, -208.185)';

	RETURN real_schema || '.' || table_name || '.' || column_name ||' datum transformed to ' || dst_srid::text;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100;
ALTER FUNCTION public.icgc_23031_to_25831(character varying, character varying, character varying)
  OWNER TO datumtest;
COMMENT ON FUNCTION public.icgc_23031_to_25831(character varying, character varying, character varying) IS 'args: schema_name, table_name, column_name - Transforms geometries in a table from EPSG:23031 to EPSG:25831 using the official catalan affine transform';

SELECT icgc_23031_to_25831('public', 'table_name', 'geom');

ST_Affine(geom, 0.9999984495712146, -0.000007587516271060413, 0.000007587516271060413, 0.9999984495712146, 129.547, 208.186)

update spatial_ref_sys set proj4text = '+proj=utm +zone=31 +ellps=intl +units=m +no_defs +nadgrids=100800401.gsb' where srid = 23031;
update spatial_ref_sys set proj4text = '+proj=longlat +ellps=intl +no_defs +nadgrids=100800401.gsb' where srid = 4230;

SELECT ST_AsEWKT(geom), ST_AsEWKT(ST_Transform(geom, 25831)) from ed50_test_points; SELECT ST_AsEWKT(geom), ST_AsEWKT(ST_Transform(geom, 23031)) from etrs89_test_points;

—
Oscar Fonts (geomati.co)