Incorrect calculation between CRS coordinate transformation

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.