Mathieu,
It is not quite true that there is no datum shift. If there was no datum
shift, you would get the following results (using EPSG:32612 to remain
in WGS84 datum, but with the UTM zone 12N similarly as EPSG:2956):
$ echo 50.937434634728845 -113.63360224940385 1005.9387180001410 |
cs2cs -d 3 EPSG:4326 EPSG:32612
314966.127 5646170.391 1005.939
The result you get in your code is the one of the below operation:
$ projinfo -s EPSG:4326 -t EPSG:2956
Candidate operations found: 2
-------------------------------------
Operation No. 1:
unknown id, Inverse of NAD83(CSRS) to WGS 84 (2) + UTM zone 12N, 1 m,
Canada - onshore and offshore - Alberta; British Columbia; Manitoba; New
Brunswick; Newfoundland and Labrador; Northwest Territories; Nova
Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon.
PROJ string:
+proj=pipeline
+step +proj=axisswap +order=2,1
+step +proj=unitconvert +xy_in=deg +xy_out=rad
+step +proj=push +v_3
+step +proj=cart +ellps=WGS84
+step +inv +proj=helmert +x=-0.991 +y=1.9072 +z=0.5129
+rx=-0.0257899075194932
+ry=-0.0096500989602704 +rz=-0.0116599432323421 +s=0
+convention=coordinate_frame
+step +inv +proj=cart +ellps=GRS80
+step +proj=pop +v_3
+step +proj=utm +zone=12 +ellps=GRS80
As the 2 CRS you use are 2D ones, the Z coordinate is left untouched
with the push/pop +v_3 steps.
If you promote the 2 CRS as 3D, you'll get
$ echo 50.937434634728845 -113.63360224940385 1005.9387180001410 |
cs2cs -d 3 EPSG:4979 "$(src/projinfo --3d EPSG:2956 -o WKT2_2019
--single-line -q)"
314967.085 5646169.595 1006.393
so you have a 3D datum shift here. Not exactly the one you expect (not
sure where your expected result comes from), but that's the best one
we've available in our database, mostly driven by what there is in EPSG.
As there are several iterations of the NAD83(CSRS) datum (the default
one used in EPSG:2956, but also NAD83(CSRS)v2, v3, v4, v5, v6, v7), it
is possible that your expected result is obtained with one of those
later datums.
API wise, you can get the 3D versions of the CRS with things like
PJ* crs_2D = proj_create(ctxt, "EPSG:2956");
PJ* crs_3D = proj_crs_promote_to_3D(ctxt, crs_2D); // in proj_experimental.h
and then use proj_create_crs_to_crs_from_pj()
Also note that the accuracy of the transformation is 1 meter, and that
anything tagged EPSG:4326 has also a potential inaccuracy of 2 meters.
Even
Le 26/11/2021 à 18:03, Mathieu Poulin a écrit :
When creating a transformation that requires a change of datum, it
does not work. For instance, going from Geographic coordinates in
WGS84 to UTM coordinates in NAD83CSRS. It does convert from geographic
lat,lon to UTM projection easting,northing but the datum remains WGS84.
I am using Proj 8.2 (but I've also tried this with Proj 7.1 and 7.2)
with C++. Here is the code with comments describing the values of the
variables, the expected output of the code and the actual output of
the code.
// ENTRY: x=-1614798.938 y=-3690228.026 z=4929942.509
GeographicLib::Geocentric::WGS84().Reverse(x, y, z, point.m_latitude,
point.m_longitude, point.m_ellipsoidal_height);
// REVERSE: point.m_latitude=50.937434634728845
point.m_longitude=-113.63360224940385
point.m_m_ellipsoidal_height=1005.9387180001410
PJ_COORD coord = proj_coord(point.m_latitude, point.m_longitude,
point.m_ellipsoidal_height, 0.0);
PJ* proj_from_WGS84_to_2956 = proj_create_crs_to_crs(0, "EPSG:4326",
"EPSG:2956", 0);
PJ_XYZ proj_point_2 = proj_trans(proj_from_WGS84_to_2956, PJ_FWD,
coord).xyz;
// EXPECTED RESULT: 314967.330, 5646169.740, 1006.379
// ACTUAL RESULT: 314967.085, 5646169.595, 1005.939
Am I doing something wrong? Am I missing something? If someone could
help me out with this bit of code, it would be hugely appreciated.
Thanks,
- Mathieu Poulin
*Mathieu Poulin*
*Développeur | Developer*
Division arpentage et géomatique
Département Recherche et Développement
Siège social
15069, boulevard Henri-Bourassa
Québec (Québec) G1G 3Z5, Canada
Tel : 1+(418) 641-0344
Sans frais: 1+ (888) 576-7898
*mvtgeosolutions.com <http://www.mvtgeosolutions.com/>*
<https://www.facebook.com/mvtgeosolutions><https://www.linkedin.com/company/mvtgeosolutions><https://www.youtube.com/channel/UCOFguYWgpbxOlKndKRqqRRQ><https://www.youtube.com/channel/UCOFguYWgpbxOlKndKRqqRRQ>
ATTENTION NOTICE
Les informations contenues dans le présent message et dans toute pièce
qui lui est jointe sont confidentielles et peuvent être protégées par
le secret professionnel, industriel et/ou droits d'auteurs. Ces
informations sont à l'usage exclusif de son ou de ses destinataires.
Si vous avez reçu ce message par mégarde, veuillez communiquer avec
l'expéditeur au +1 (888) 576-7898 poste 101, l'effacer de tout disque
dur ou autre média sur lequel il peut être enregistré et ne pas en
conserver de copie. Merci.
This E-mail message and any attachment thereto contain confidential
information which may be privileged, industrial secrets or copyrighted
and which is intended for the exclusive use of its addressee(s). If
you have received this communication in error, please immediately
notify us by telephone at +1 (888) 576-7898 poste 101, erase it from
any hard disk or other medium on which it may have been saved and do
not keep any copy thereof. Thank you.
_______________________________________________
PROJ mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/proj
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
PROJ mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/proj
