¿Conversión de un sistema de coordenadas personalizado oblicuo y no cartesiano con PostGIS?

Question

Tengo un sistema de coordenadas personalizado que es:

• No cartesiano : La unidad del eje X es 25m Eje Y 6.25m
• Oblicuo : Ninguno de los azimuts de ambos ejes es Norte-Sur ( 0° ). El eje X va 54.2° Eje Y 324° .
• Ortogonal : Allí 90° entre ambos ejes.

Por ejemplo si tengo esos 3 puntos (La primera columna es long/lat, la segunda está en el sistema de coordenadas personalizado):

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)

¿Puedo escribir un PROJ.4 texto para utilizar con ST_Transform ?

Inspirado en el Ejemplo de documentación PostGIS ¿funcionaría algo así?

 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)

Vi esto buena respuesta para la parte oblicua pero es con una rejilla métrica cartesiana.

Edita:

He simplificado para el ejemplo, pero la referencia espacial real es:

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)

La esquina más oriental es POINT(901 501) . Está situada en el Caribe, entre México y Cuba.

Ya sé cómo convertirlo en mi aplicación pero me gustaría hacerlo en PostGIS almacenando toda la referencia espacial en spatial_ref_sys.proj4text y utilizando ST_Transform. Podría ser una transformación de la rejilla personalizada a SRID 4326 o de la cuadrícula personalizada a coordenadas geográficas.

Answer 1

Para su problema en particular, parece que lo mejor es utilizar ST_Affine() para convertir desde el sistema local a EPSG:3795, a partir del cual podrá utilizar ST_Transform() para llevarlo a cualquier otro sistema. El problema es que no parece haber ninguna herramienta fácil para convertir un conjunto de tres pares de puntos en una transformación afín. ¡Qué locura! Tuve que escribir un script en python para hacerlo, y generé la siguiente transformación para tus pares de puntos (usando los tres primeros, y probado aquí con el cuarto):

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)

Estos parámetros afines funcionarán para cualquier punto de su malla local y los convertirán a coordenadas EPSG:3795. En combinación con ST_Transform() puedes dar un paso más hacia la geografía:

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)

Creo que necesito convertir mi programa python en un script PL/PgSQL para que otros puedan hacer este truco un poco más fácil.

Answer 2

Proyecto.4

No es posible 1 , 2 para hacer una transformación afín con proyecto.4 porque sólo hay un parámetro de escala para todas las dimensiones ( +k_0 ) y +towgs84 es sólo un transformación de siete parámetros (3 traslaciones + 3 rotaciones + 1 escalado).

Ha sido propuesta comme +xform pero aún no se ha aplicado.

Texto conocido

No es posible con WKT versión 1 .

El nuevo texto conocido versión 2 (WKT 2 - ISO 19162:2015) lo hace posible introduciendo descripciones de ingeniería derivada CRS que puede tener sistemas de coordenadas afines .

La conversión a proyecto.4 es aún no se ha aplicado .

PostGIS

En lugar de almacenar un proyecto.4 text cadena en la columna spatial_ref_sys.proj4text Es posible almacenar la transformación en otro lugar como text que describe la matriz de transformación y utilizar ST_Affine(geom, a, b, d, e, xoff, yoff) para ejecutarlo.

Parámetros de transformación

Ecuación matricial para 3 puntos A, B y 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│

Se puede resolver calculando │P│ -1 .│P'│ = │T│

Calcular la matriz de transformación en 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]]

Puede ejecutarlo en línea y modifícalo con tus propios parámetros.

Solución

│ 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)

El resultado coincide con el 4º punto de la pregunta.

Almacenar los parámetros de transformación en una tabla

Así que puede utilizarlo de forma similar a 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)

Convertir a cualquier otro sistema de coordenadas

Por ejemplo, a coordenadas geográficas de longitud/latitud utilizando ST_Transformación :

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

Esta pregunta es casi la misma que " Utilización de la biblioteca PROJ.4 para transformar coordenadas de un sistema de coordenadas local a un sistema de coordenadas global utilizando puntos de control terrestres. ", pero utilizando PostGIS.

Los puntos pueden transformarse utilizando una matriz de transformación afín 2D mediante ST_Affine que tiene seis coeficientes, que puedes obtener resolviendo con Matlab, o mostrando a continuación con 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

Copie la salida de los seis coeficientes para usarla con SQL, así:

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)

Tenga en cuenta que ST_SnapToGrid se utilizó para redondear los pequeños errores de coma flotante. Además, he desajustado el SRID (a 0), ya que las coordenadas ya no son EPSG:3795.

No creo que PROJ.4 pueda definir proyecciones hacia delante/inversas basadas en coeficientes de transformación afín.