[PROJ] Change of datum not working
Even Rouault
even.rouault at spatialys.com
Fri Nov 26 09:32:06 PST 2021
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
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj
--
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20211126/e8c4922b/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Outlook-31dfdvgz.png
Type: image/png
Size: 8696 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20211126/e8c4922b/attachment-0004.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Outlook-gxopdjhy.png
Type: image/png
Size: 1444 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20211126/e8c4922b/attachment-0005.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Outlook-aio4tfz2.png
Type: image/png
Size: 880 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20211126/e8c4922b/attachment-0006.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Outlook-ccv4tpf3.png
Type: image/png
Size: 1584 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20211126/e8c4922b/attachment-0007.png>
More information about the PROJ
mailing list