Converting from oblique and non-cartesian custom coordinates system with PostGIS?

Question

I have a custom coordinates system that is:

• Not cartesian: The X-axis unit is 25m & Y-axis 6.25m
• Oblique: None of both axis' azimuth is North-South (0°). X-axis goes 54.2° & Y-axis 324°.
• Orthogonal: There 90° between both axis.

For instance if I have those 3 points (The first column is long/lat, the second is in the custom coordinates system):

SRID=4326;POINT(-85.647 21.318) <=> POINT(0 0)
SRID=4326;POINT(-85.599 21.353) <=> POINT(0 1000)
SRID=4326;POINT(-85.794 21.497) <=> POINT(1000 0)

Can I write a custom PROJ.4 text to use with ST_Transform?

Inspired from the PostGIS documentation's example, would something like that work?:

 SELECT ST_AsText(
   ST_Transform(
     ST_GeomFromText('POINT(1000 1000)'), custom_proj4, 4326))
 FROM (
   SELECT
     '?????????????????????????????????'::text AS custom_proj4
 ) AS data;
                                  st_astext


POINT(-85.7467905609375 21.5321592956322)

I saw this good answer for the oblique part but it is with a cartesian metric grid.

Edit:

I simplified for the example, but the real spatial reference is:

SRID=3795;POINT( 7217.670 193658.490) <=> POINT( 901   501)
SRID=3795;POINT(-7052.270 213450.630) <=> POINT(1877   501)
SRID=3795;POINT(66458.500 266451.180) <=> POINT(1877 15001)
SRID=3795;POINT(80728.440 246659.040) <=> POINT( 901 15001)

The most eastern corner is POINT(901 501). It is located in the caribbean, between Mexico and Cuba.

I already know how to convert it in my application but I'd like to do it in PostGIS by storing the whole spatial reference in spatial_ref_sys.proj4text and using ST_Transform. It could be a transformation from the custom grid to SRID 4326 or from the custom grid to geographical coordinates.

Answer 1

For your particular problem, it seems like the nicest thing to do is to use ST_Affine() to convert from the local system into EPSG:3795, from which you can then use ST_Transform() to take it to any other system. The trouble is, there seems to be no easy tooling to convert a set of three point-pairs into an affine transformation! Crazy! I had to write a python script to do it, and generated the following transform for your point pairs (using the first three, and tested here on the fourth):

select st_astext(st_affine(
  'POINT(901 15001)'::geometry, 
  -14.620840163934428, 5.0697082758620695, 
  20.278831967213137, 0.6552103448275854, 
  17851.12314149802, 173556.00201478234));

         st_astext         
---------------------------
 POINT(80728.44 201656.04)
 (1 row)

These affine parameters will work for any point in your local grid and convert them to EPSG:3795 coordinates. In combination with ST_Transform() you can one-step out to geographics:

select st_astext(st_transform(st_setsrid(st_affine(
  'POINT(901 15001)'::geometry, 
  -14.620840163934428, 5.0697082758620695, 
  20.278831967213137, 0.6552103448275854, 
  17851.12314149802, 173556.00201478234), 3795), 4326));

         st_astext         
---------------------------
 POINT(-85.0487027252303 21.589299875787)
 (1 row)

I think I need to convert my python program into a PL/PgSQL script so others can do this trick a little easier.

Answer 2

proj.4

It is not possible ^1,2 to do an affine transformation with proj.4 because there is only one scaling parameter for all dimensions (+k_0) and +towgs84 is only a seven-parameter transformation (3 translations + 3 rotations + 1 scaling).

It has been proposed as +xform but not implemented yet.

Well-known text

It is not possible with WKT version 1.

The new Well-known text version 2 (WKT 2 - ISO 19162:2015) makes it possible by introducing descriptions of derived engineering CRS that can have affine coordinate systems.

The conversion to proj.4 is not implemented yet.

PostGIS

Instead of storing a proj.4 text string in the column spatial_ref_sys.proj4text It is possible to store the transformation somewhere else as text describing the transformation matrix and use ST_Affine(geom, a, b, d, e, xoff, yoff) to execute it.

Transformation parameters

Matrix equation for 3 points A, B & C:

│P│.│T│ = │P'│

│Ax Ay 1│ │   a    d 0│   │Ax' Ay' 1│
│Bx By 1│.│   b    e 0│ = │Bx' By' 1│
│Cx Cy 1│ │xoff yoff 1│   │Cx' Cy' 1│

It can be solved by calculating │P│^-1.│P'│ = │T│

Calculate the transformation matrix in Ruby

│ORIGINAL│.│TRANSFORMATION│ = │TRANSFORMED│

require "matrix"

ORIGINAL = Matrix[
                   [ 901,   501, 1],
                   [1877,   501, 1],
                   [1877, 15001, 1]
                 ]

TRANSFORMED = Matrix[
                      [ 7217.67, 193658.49, 1],
                      [-7052.27, 213450.63, 1],
                      [66458.50, 266451.18, 1]
                    ]

puts TRANSFORMATION = ORIGINAL.inverse * TRANSFORMED

=> Matrix[[-14.620840163934428, 20.278831967213137, 0/1], [5.0697082758620695, 3.6552103448275854, 0/1], [17851.12314149802, 173556.00201478234, 1/1]]

You can run it online and modify it with your own parameters.

Solution

│ 901   501 1│ │  -14.62     20.28 0│   │ 7217.67 193658.49 1│
│1877   501 1│.│    5.07      3.66 0│ = │-7052.27 213450.63 1│
│ 901 15001 1│ │17851.12 173556.00 1│   │80728.44 246659.04 1│

SELECT ST_AsEWKT(
         ST_SetSRID(
           ST_Affine(geom, t[1], t[2], t[3], t[4], t[5], t[6]), 
         to_srid)
       )
FROM (SELECT 
        'POINT(1877 15001)'::geometry AS geom,
        3795 AS to_srid,
        array[-14.620840163934428, 5.0697082758620695, 
               20.278831967213137, 3.6552103448275854, 
                17851.12314149802, 173556.00201478234] AS t
      ) AS parameters;

                        st_asewkt
----------------------------------------------------------
SRID=3795;POINT(66458.5 266451.18)

The result matches the 4th point of the question.

Store the transformation parameters in a table

So you can use it in a way similar to spatial_ref_sys.

CREATE TABLE affine_transformations
(
  id integer,
  to_srid integer,
  t numeric[6]
);

INSERT INTO affine_transformations 
VALUES (1, 3795, '{-14.620840163934428, 5.0697082758620695, 
                    20.278831967213137, 3.6552103448275854, 
                     17851.12314149802, 173556.00201478234}');

SELECT ST_AsEWKT(
         ST_SetSRID(
           ST_Affine('POINT(1877 15001)', t[1], t[2], t[3], t[4], t[5], t[6]), 
         to_srid)
       )
FROM affine_transformations
WHERE id = 1;

                        st_asewkt
----------------------------------------------------------
SRID=3795;POINT(66458.5 266451.18)

Convert to any other coordinate system

For example to geographic longitude/latitude coordinates using ST_Transform:

SELECT ST_AsEWKT(
         ST_Transform(
           ST_SetSRID(
             ST_Affine('POINT(1877 15001)', t[1], t[2], t[3], t[4], t[5], t[6]), 
           to_srid),
         4326)
       )
FROM affine_transformations
WHERE id = 1;

                        st_asewkt
----------------------------------------------------------
SRID=4326;POINT(-85.2039368976408 22.1707571653241)

Answer 3

This question is nearly the same as "Using PROJ.4 library to transform from local coordinate system coordinates to global coordinate system using ground control points?", but using PostGIS.

The points can be transformed using a 2D affine transformation matrix via ST_Affine, which has six coefficients, which you can get by solving using Matlab, or shown below with Numpy (Python):

import numpy as np

# input and output coordinates, with 1 padded at the end
X = np.array([
    (7217.67, 193658.49, 1),
    (-7052.27, 213450.63, 1),
    (66458.5, 266451.18, 1),
    (80728.44, 246659.04, 1),
    ], dtype='d')
Y = np.array([
    (901, 501, 1),
    (1877, 501, 1),
    (1877, 15001, 1),
    (901, 15001, 1),
    ], dtype='d')

# Make sure we have the same number in each array
assert X.shape == Y.shape

# Use a least-squares solution to find coefficients for the linear system of equations
A, res, rank, s = np.linalg.lstsq(X, Y)

# Make sure this is a small number
print('Maximum residule was: ' + str(np.abs(res).max()))

# Get the coefficients for PostGIS' ST_Affine
a, d, _, b, e, _, xoff, yoff, _ = A.flatten()
print(', '.join(repr(x) for x in [a, b, d, e, xoff, yoff]))
# -0.023393344939273676, 0.032446131207436515, 0.12978451755815532, 0.093573371978935407, -5213.6233320082156, -18557.019740500145

Copy the output of the six coefficients for use with SQL, like this:

SELECT ST_AsEWKT(g) || ' <=> ' || ST_AsEWKT(
  ST_SnapToGrid(ST_SetSRID(ST_Affine(g,
      -0.023393344939273676, 0.032446131207436515,
       0.12978451755815532, 0.093573371978935407,
      -5213.6233320082156, -18557.019740500145),
  0), 1e-7)) AS result
FROM (
  SELECT 1 AS id, 'SRID=3795;POINT( 7217.670 193658.490)'::geometry AS g
  UNION SELECT 2, 'SRID=3795;POINT(-7052.270 213450.630)'
  UNION SELECT 3, 'SRID=3795;POINT(66458.500 266451.180)'
  UNION SELECT 4, 'SRID=3795;POINT(80728.440 246659.040)'
) f
ORDER BY id;
                          result
----------------------------------------------------------
 SRID=3795;POINT(7217.67 193658.49) <=> POINT(901 501)
 SRID=3795;POINT(-7052.27 213450.63) <=> POINT(1877 501)
 SRID=3795;POINT(66458.5 266451.18) <=> POINT(1877 15001)
 SRID=3795;POINT(80728.44 246659.04) <=> POINT(901 15001)
(4 rows)

Note that ST_SnapToGrid was used to round off small floating point errors. Also, I've unset the SRID (to 0), since the coordinates are no longer EPSG:3795.

I don't think PROJ.4 can define forwards/inverse projections based on affine transformation coefficients.