5

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 (). 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.

Graphical representation of the

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.

  • 1
    If you combine the oblique answer with an ST_Scale to push your coordinates into a metric cartesian space first, you could be OK. – Paul Ramsey Jan 25 '17 at 18:24
  • @PaulRamsey That's a good idea, sure it should work. But then I still need something out of PostGIS to calculate the scaling parameters. And doing the transformation in 3 steps (calculate scaling, ST_Scale, ST_Transform) sounds like a workaround. I'm quite sure it's possible using only PROJ.4 but the documentation is quite poor (with stuff "grepped out" from directories). – Victor Jan 25 '17 at 18:48
  • Are you sure about your corner ticks? The "positive X" one is actually to the west of the origin... which would imply a southerly oriented map, except that the "positive Y" example is north of the origin, as one might expect. – Paul Ramsey Jan 25 '17 at 20:08
  • They also appear to be in the Caribbean Sea – Paul Ramsey Jan 25 '17 at 20:15
  • @PaulRamsey Edited with the real coordinates. It is indeed in the Caribbean Sea. It an ad hoc local grid for a specific research survey. – Victor Jan 26 '17 at 18:28
2
+50

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.

2

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.

  • Thanks. In fact I'm woking on exotic SRID 3795 (neither UTM or long/lat) and converted it to 4326 to make a general example. So I'll check my calculations, starting with the azimuth. But dont you think it's possible to scale with PROJ.4? +proj=omerc has to be cartesian? – Victor Jan 25 '17 at 19:06
  • Scaling is possible, but not independent of axis as far as I can tell: proj4.org/parameters.html. If both axes are the same the +to_meter parameter can be used to scale them into meters. – Paul Ramsey Jan 25 '17 at 20:01
  • I did geometric calculations, not geographic. It was a bad idea of mine to convert to long/lat for the example. That probably explains the difference you find in azimuth. I suppose the difference in the positions you find are due to the fact your PROJ.4 origin is geographic: +lat_0=21.318 +lonc=-85.647. – Victor Jan 26 '17 at 18:35
  • Nope, there seems to be no way to add an affine transformation step to a proj4 string to achieve what you want. But since PostGIS supports affines, you can still do it all in PostGIS. Why the obsession with a pure-proj solution? Simplest solution is to affine from your local grid into ESPG:3795, then just proj4 into whatever output projection you desire. – Paul Ramsey Jan 26 '17 at 18:46
  • I thought that if I could store the whole spatial reference in spatial_ref_sys.proj4text, I could then shove raw coordinates into PostGIS and it would take care of everything. If I do the scaling in a different step, I sill need the application to calculate the scale and store it in some other place than spatial_ref_sys. – Victor Jan 26 '17 at 18:55
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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.