Incorrect calculation between CRS coordinate transformation

Hello,

I want to convert the point 4142408.339 mN, 368674.783mE (Zone 19F or 19S) in EPSG:2084 to a WGS-84 system. I choose EPSG:32719 as it is WGS-84 for 19S. I ahve this script that correctly transform to lat long, but incorrectly calculate the coordinates in the new CRS:

% Define point in EPSG 2084
y1 = 4142408.339;
x1 = 368674.783;
% Define the initial (1) and final (2) CRS
p1 = projcrs(2084); % "Hito XVIII 1963 / UTM zone 19S"
p2 = projcrs(32719); % "WGS 84 / UTM zone 19S"
% Transform the coordinates from EPSG 2084 to EPSG 32719
[lat, lon] = projinv(p1, x1, y1); % Convert to geographic coordinates (latitude and longitude)
[x2, y2] = projfwd(p2, lat, lon); % Convert from geographic coordinates to final CRS
% Display the results
fprintf('Coordinates with EPSG:32719: (%f, %f)\n', y2, x2);

Coordinates with EPSG:32719: (4142531.755753, 368681.136119)

From GPS confirmation and topography, we know the point is 4142534.593 mN and 368760.170 mE. However, MATLAB answer is 4142531.755753 mN and 368681.136119 mE.

What would bethe error?

Thank you!

Answer 1

[Umar]

Hi @Cristobal,

I reviewed the workflow for converting coordinates from EPSG:2084 (Hito XVIII 1963 / UTM 19S) to EPSG:32719 (WGS84 / UTM 19S). The difference you observed (~70–80 m in easting) is expected — the source of the discrepancy is a datum mismatch, not a coding mistake.

Explanation:

• `projinv(p1, x, y)` returns latitude/longitude in the geographic CRS associated with EPSG:2084 — i.e., Hito XVIII 1963, not WGS84 ([projinv docs]( https://www.mathworks.com/help/map/ref/projcrs.projinv.html )).

• `projfwd(p2, lat, lon)` then assumes those coordinates are already in WGS84, since EPSG:32719 is tied to WGS84. MathWorks explicitly warns that if the geographic CRS of the input does not match the target CRS, the result may be inaccurate ([projfwd docs]( https://www.mathworks.com/help/map/ref/projcrs.projfwd.html )).

• MATLAB’s Mapping Toolbox performs projection math but does *not automatically apply datum transformations* ([example: transforming CRSs]( https://www.mathworks.com/help/map/transform-coordinates-to-different-pcrs.html )).

This explains why the MATLAB result differs from GPS/topographic measurements. The missing step is applying the datum transformation from Hito XVIII 1963 → WGS84.

Options to resolve:

1. Inside MATLAB: Apply the EPSG-defined transformation parameters manually, then project to EPSG:32719. 2. Outside MATLAB: Use PROJ/GDAL or QGIS, which apply datum transformations automatically, e.g.:

   cs2cs +init=epsg:2084 +to +init=epsg:32719

Appendix: EPSG Transformation Parameters (EPSG:1529)

• dX = +18.38 m
• dY = +192.45 m
• dZ = +96.82 m
• rotX = +0.056 arc‑seconds
• rotY = –0.142 arc‑seconds
• rotZ = –0.200 arc‑seconds
• scale = –0.0013 ppm

(Source: [EPSG:1529]( https://epsg.io/1529 ))

MATLAB Snippet

function [x_wgs, y_wgs] = convert_HitoXVIII_to_EPSG32719(x_hito, y_hito)
  p1 = projcrs(2084); 
  [lat_hito, lon_hito] = projinv(p1, x_hito, y_hito);

    [Xh, Yh, Zh] = geodetic2ecef(referenceEllipsoid('International1924'), ...
                                 lat_hito, lon_hito, 0);

    dX=18.38; dY=192.45; dZ=96.82;
    rx=0.056; ry=-0.142; rz=-0.200; % arc-seconds
    s=-0.0013; % ppm
    arcsec2rad = pi/(180*3600);
    R = eye(3) + [0 -rz*arcsec2rad ry*arcsec2rad;
                  rz*arcsec2rad 0 -rx*arcsec2rad;
                 -ry*arcsec2rad rx*arcsec2rad 0];
    scale = 1 + s*1e-6;
    vec = scale * (R * [Xh;Yh;Zh]) + [dX;dY;dZ];
    [lat_wgs, lon_wgs, ~] = ecef2geodetic(referenceEllipsoid('WGS84'), ...
                                          vec(1), vec(2), vec(3));

    p2 = projcrs(32719); 
    [x_wgs, y_wgs] = projfwd(p2, lat_wgs, lon_wgs);
  end

Applying this procedure should produce coordinates consistent with GPS/topographic measurements.

[Torsten]

Only closer ...

% Define point in EPSG 2084
x_hito = 368674.783;
y_hito = 4142408.339;
format long
[x_wgs, y_wgs] = convert_HitoXVIII_to_EPSG32719(x_hito,y_hito)

x_wgs = 
     3.687611121856684e+05

y_wgs = 
     4.142540017764304e+06

function [x_wgs, y_wgs] = convert_HitoXVIII_to_EPSG32719(x_hito, y_hito)
  p1 = projcrs(2084); 
  [lat_hito, lon_hito] = projinv(p1, x_hito, y_hito);
  [Xh, Yh, Zh] = geodetic2ecef(referenceEllipsoid('International1924'), ...
                               lat_hito, lon_hito, 0);
  dX=18.38; dY=192.45; dZ=96.82;
  rx=0.056; ry=-0.142; rz=-0.200; % arc-seconds
  s=-0.0013; % ppm
  arcsec2rad = pi/(180*3600);
  R = eye(3) + [0 -rz*arcsec2rad ry*arcsec2rad;
                rz*arcsec2rad 0 -rx*arcsec2rad;
               -ry*arcsec2rad rx*arcsec2rad 0];
  scale = 1 + s*1e-6;
  vec = scale * (R * [Xh;Yh;Zh]) + [dX;dY;dZ];
  [lat_wgs, lon_wgs, ~] = ecef2geodetic(referenceEllipsoid('WGS84'), ...
                                        vec(1), vec(2), vec(3));
  p2 = projcrs(32719); 
  [x_wgs, y_wgs] = projfwd(p2, lat_wgs, lon_wgs);
end

[Cristóbal]

Thank you @Umar, I made a reply below

Answer 2

[Cristóbal]

Hello Umar,

Your explanation of the issue was very clear.

The script works perfectly. I did a slightly modification, as EPSG defines three types of transformations with different values of dX, dY, dZ and rx,ry,rz and s for going from "Hito XVII 1963" to "WGS 84".

The one used to get the values is "Chile - Tierra del Fuego onshore. , accuracy 44.0 m, code 1892 (default) [3]" which have the values dX=16, dY=196, dZ=93 and rx=0, ry=0, rz=0, s=0.

https://epsg.io/2084-1892

https://epsg.io/1892

With that, I get the proper coordinates y=4142534.593804 and x=368760.169510 which are almost exact (difference less than 1 mm). Difference may be due to decimals in the original coordinates 52°51'3.3064"S / 70°57'0.5434"W, transformed to UTM: 4142408.339 mN and 368674.783 mE

If you agree with me, I will accept this answer. Thank you for your answer!

[Cristóbal]

Umar,

Is it possible to use the script for several values? I've got this error:

Error using *

Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element of the matrix

individually, use TIMES (.*) for elementwise multiplication.

Thank you!

[Umar]

Hi @Cristóbal,

I'm glad the solution worked well for you! However, I want to make sure credit goes where it's due - @Torsten deserves the primary recognition here. While I provided the theoretical explanation and initial code framework, @Torsten was the one who actually implemented, tested, and validated the solution with real data. His mathematical expertise and practical validation were crucial to making this work.

Regarding your question about using the script for multiple values - that error suggests you're trying to pass arrays/matrices where the function expects scalar values. The current function is designed for single coordinate pairs. To handle multiple points, you'd need to modify it to loop through the coordinates or vectorize the operations properly.

Thanks for acknowledging the collaborative effort - that's what makes these mathworks technical communities so valuable.

[Cristóbal]

Thank you very much for your support @Torsten... It provided a valuable solution that helped me a lot, as well as a brilliant code. Thank you again @Umar, I was able to create a for loop to do the calculations for a batch of values. I will proceed to accept the solution.