¿Cómo puedo determinar si una transformación entre dos SIR es simplemente una escala/traducción o algo más complejo?

Question

Estamos utilizando pyproj para transformar entre sistemas de referencia de coordenadas. ¿Cuál es una buena manera de comprobar si la transformación es simplemente escala/translación/unidades? Básicamente, quiero saber si una cuadrícula en el primer crs seguirá siendo en una cuadrícula en el crs de destino.

Edición: Quiero añadir algo de claridad sobre lo que estoy tratando de detectar. Supongamos que tengo una cuadrícula de 3x4 de 12 puntos representados por 3 x y 4 y valores. Cuando transformo estos puntos, a veces la cuadrícula se distorsiona de tal manera que 12 x,y se requieren pares. Pero a veces, la transformación da lugar a una cuadrícula de 3x4 que todavía puede representarse utilizando 3 x y 4 y valores. Son estos últimos casos los que queremos detectar.

Un ejemplo es EPSG:4236 a EPSG:3395 . Otras transformaciones distorsionan la red, como EPSG:4236 a EPSG:2193 .

He enumerado un montón de detalles a continuación, pero creo que esta va a ser una pregunta bastante fácil. Lo ideal es que haya algo sencillo en el pyproj.Transformer que pueda inspeccionar, o una forma fácil de comparar las definiciones de las proyecciones en algún formato.

1. Prueba empírica

   Obviamente, puedo transformar algunos puntos y luego comprobar si siguen siendo cuadriculados. Me gustaría encontrar una forma mejor.

2. Comparar DATUM

   Utilizando pyproj.CRS('EPSG:4236').datum los datos son

   • EPSG:4236 : DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]], ID["EPSG",6326]]
   • EPSG:3395 : DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]], ID["EPSG",6326]]
   • EPSG:2193 : DATUM["New Zealand Geodetic Datum 2000", ELLIPSOID["GRS 1980",6378137,298.257222101, LENGTHUNIT["metre",1]], ID["EPSG",6167]]

   Así que los DATUM de 4326 y 3395 son exactamente los mismos. ¿Es el caso que dos crs con la misma cadena de datos exacta es siempre una escala/traducción? ¿Es suficiente con comparar sólo la parte del DATUM ( World Geodetic System 1984 vs New Zealand Geodetic Datum 2000 ) o sólo la parte de la identificación ( EPSG 6326 vs EPSG 6167 )? O bien, ¿esto no va a funcionar en absoluto (pueden dos SIR compartir el mismo DATUM pero aún así requerir la distorsión al transformarse)?

   Además, aunque eso funcione, ¿tiene cada crs un "dato"? Sería bueno detectar esto incluso con crs personalizados.

3. Comparar PROJ4

   Las cadenas del proyecto4:

   • EPSG:4326 : +proj=longlat +datum=WGS84 +no_defs +type=crs
   • EPSG:3395 : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs
   • EPSG:2193 : +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs +type=crs

   ¿Puedo hacer algo con esto? Parece que no hay suficiente información, aquí, pero tal vez algunos valores tienen "valores por defecto"? ¿Qué partes de estas cadenas indican el escalado/transformación en contraposición a otros cálculos?

4. Comparar WKT

   Las cadenas WKT o las representaciones JSON tienen los mismos parámetros básicos, pero son más detalladas y pueden ser más sencillas de utilizar. No las publicaré aquí del todo. Pero, por ejemplo, EPSG:3395 tiene EPSG:4326 como base_crs pero, por supuesto, EPSG:2193 no lo hace. ¿Compararían base_crs ¿es suficiente?

5. Inspeccionar el pyproj Transformer objeto

   Estoy usando un pyproj Transformer Así que tal vez lo mejor sea simplemente inspeccionarla.

   • Uso del transformador definition directamente podría ser mejor que intentar comparar las cadenas de proj4; ¿qué parte(s) debo utilizar?

   • Como alternativa, el operations está vacía para la transformación de 4236 a 3395, pero no para la de 4326 a 2193. ¿Es tan sencillo como eso?

     t = pyproj.Transformer.from_crs('EPSG:4236', 'EPSG:3395') t.definition 'proj=pipeline step proj=axisswap order=2,1 step proj=unitconvert xy_in=deg xy_out=rad step proj=merc lon_0=0 k=1 x_0=0 y_0=0 ellps=WGS84' t.operations ()

     t = pyproj.Transformer.from_crs('EPSG:4236', 'EPSG:2193') t.definition 'proj=pipeline step proj=axisswap order=2,1 step proj=unitconvert xy_in=deg xy_out=rad step proj=tmerc lat_0=0 lon_0=173 k=0.9996 x_0=1600000 y_0=10000000 ellps=GRS80 step proj=axisswap order=2,1' t.operations (<Coordinate Operation: Transformation> Inverse of NZGD2000 to WGS 84 (1) Area of Use:

     • name: New Zealand
     • bounds: (160.6, -55.95, -171.2, -25.88), <Coordinate Operation: Conversion> New Zealand Transverse Mercator 2000 Area of Use:
     • name: New Zealand - onshore
     • bounds: (166.37, -47.33, 178.63, -34.1))

Answer 1

Como respuesta rápida, la mayoría de las transformaciones son más complejas que una escala/translación/unidad, por lo que se buscan excepciones. La distorsión proviene del cambio de dato (es decir, la transformación, de lat/long a lat/long) y del tipo de proyección (de lat/long a XY).

Tener el mismo (o girado alrededor de la vertical) punto de referencia es un requisito previo, pero si tiene diferentes proyecciones sobre el punto de referencia, esto "distorsionará" su cuadrícula regular.

Edición: como observación, tus cuadrados en EPSG4326 se convertirán en rectángulos en epsg3395 (alargamiento dependiendo de la latitud). Si no te importan esos diferentes reescalados en las direcciones X e Y, entonces todos los proyecciones cilíndricas con un eje vertical dará lugar a una retícula de rectángulos (cuadrados en el caso especial del Plate Carrée). Lamentablemente, la denominación "cilíndrica" no siempre figura en la descripción del SIR, por lo que habría que construir primero una tabla de búsqueda).