A set of functions to calculate coordinate transformations between different reference ellipsoids and different projections, including tools on:
- transformation from cartesian to geographic coordinates and back
- transformation from geographic coordinates to transverse mercator mapping or lambert conformal conical mapping and back
- transformation from geographic to UTM and back which can handle irregular zones and pole mapping
- 3D/2D/1D similarity transformation (Helmert transformation)
- determination of the parameters of a 3D/2D/1D-Helmert transformation
- applying residual corrections after performing a helmert transformation
- reading and using NTv2 transformation parameters
- 3D transformation between ITRS and ETRS frames
- 3D/2D affine transformation and its parameter determination
- 3D/2D to 2D projective transformation and its parameter determination
- Molodensky transformation
Some projections and ellipsoids are already defined in mat-Files.
This functions are often used in geodesy when transforming between different coordinate systems, e.g. UTM to GK or GPS-data (WGS84) to UTM.
Although some might be incorporated in the mapping toolbox, at least I don't have that...
An example document with background information is added.
A great acknowledgement to Andrea Pinna of the University of Cagliari who improved the helmert parameter estimation functions in speed and memory usage!
Please notify me of any bugs you may encounter.
Peter Wasmeier (2020). Geodetic Transformations (https://www.mathworks.com/matlabcentral/fileexchange/9696-geodetic-transformations), MATLAB Central File Exchange. Retrieved .
Many thanks, Peter!
the lambertcc2ell function called with the 'france_2_et' projection parameter cannot transform NTF to WGS84 coordinates as the lambert projection for NTF is defined on Clarke 1880 ellipsoid rather than on WGS84 ellipsoid.
A direct unmapping to ellipsoidal WGS84 coordinates is not possible, but French cadastral service provides a simple translation for the ellipsoid change between Clarke and WGS84: [Tx, Ty, Tz] = [-168 -60 320]
Then, your workflow would be like
k1=lambertcc2ell([600848.20, 2428281.98], 'france_2_et') % Unmapping to ellipsoidal Clarke1880 coordinates
k2=ell2cart(k1,'clarke1880') % Transform ellipsodal to cartesian Clarke1880
k3=k2+[-168 -60 320] % Change ellipsoid from Clarke1880 to WGS84
k4=cart2ell(k3,'wgs84') % Transform cartesian to ellipsoidal WGS84
k4 = [2.34805736708604 48.8528806028417 43.157798551023]
Thanks for the quick reply, Peter.
I may be using lambertcc2ell incorrectly. Here's an example in Paris:
>> lambertcc2ell([600848.20, 2428281.98], 'france_2_et')
ans = 2.3487804, 48.8529498
But we obtain approximately
when we use the same "Lambert II Etendu to WGS84" conversion in either of these websites:
The discrepancy between these websites is of the order of one meter; while it is about 51 meters for lambertcc2ell in the X direction (longitude) and 7 meters in the Y direction (latitude).
the lambertcc2ell function has been tested with some sample data from France which showede to be calcualted correctly.
Please let me know about the websites and the samples used there to control the calculations.
Has lambertcc2ell([x ; y], 'france_2_et') been sufficiently tested?
Comparing its output to outputs from some websites (e.g. in the Paris area), the function seems to have some inaccuracy.
Hi Peter ! Thank you for your contribution. I am novice in Datum and Transformation, and I am facing a problem of transforming a waypoints(lat,long) from NAD27 to WGS84. I know the Zone and the Ellipsoid coefficient. I read your pdf, but I am lost in all the details.Where do I get all the parameters you are using in your example? I need +/-1 meter accuracy. Thank you for your help.
Really nice work!!!
you are right, there was a typo in the derivative you mentioned.I corrected that.
The treason why both coordinate sets get reduced by the same centroid coordinates is due to the 10 parameter transformation. If there were two sets of reduction translations, it would become a 13 parameter transformation otherwise. Generally, the shift parameters do not necessarily need to be the centroid of any of the coordinate sets...
Hello Dr. Wasmeier,
your wellwritten toolbox helped me alot! As a geodesy student from the "Hochschule Neubrandenburg" I used your helmert3d.m as an ispiration for a conditioned Gauß-Helmert-transformation ,wich I need for my Master-thesis. I found two suspicios things during my programming:
in Line 217, wich is A(1:3:end,6) = dF_X / d_kappa:
in Col 140 [the differential m*(z-zs)*r13], you accidentally put cosin_omega instead of cosin_kappa, wich results in slightly wrong kappas for the first iteration of your algorithm. Those "wrong" Kappas will use some "room" within the next least-square-iteration. Whatsoever... With my Example-Data this error didn't do any significant Damage :)
The second thing is that I don't really understand the build of your model:
if it is build like (line 230):
it is an efficient move to calculate the 'short-observations' right away, but it looks like you are reducing both: the origin and the destination coordinates with the same core-reduction values...
Obviously it is working... and it might just be my lack of imagination, but I chose to reduce the origin_coordinates with the origin-core-reduction and the destination_coordinates with the destination-core-reduction. Whatsoever... my calculation deviaties in arround e-8 [m] wich might be due to the "kappa Error".
Thanks again for this awesome toolbox! If there is any chance that you have written a code for a 3D Helmert-transformation with fixed omega- and phi-angles.... let me now. It is a great help to verrify my work :)
really thankful for this good code, and hope to write another in geodesy.
Fantastic! Well-written, with excellent tutorial documentation.
Michael, of course you are right. I corrected the wrong file. Thanks!
There seems to be a small bug in d2affinetrafo.m: In line 46f. you check for the parameter vector length being equal to 6. But in lines 68 and 70 you adress p(6:7). Invaluable toolbox nonetheless!
Very useful. Thank you.
Concerning the definition of the coordinate system "helmert3D" is based on:
I suppose it is a right-handed system with the directions of the rotations defined as following:
x-axis: from z to y
y-axis: from x to z
z-axis: from y to x
Is that correct?
Genau das habe ich gesucht. Grossartig!
i think it is good
Great tools, but the projection.mat file seems to be corrupted - or is that just my download?
for help me
In itrstrafo function, ITRF2014 was added as well as transformations frim IGS2008 to NAD83.
There was a bug when one tried to enter parameters from the pre-defined parameter list (not the number of parameters was chekced, but the number of chjaracters of the parameter name).
I corrected a minor typo in helmert3d function and added a short comment in the manual regarding the usage of PROJ.4 transformation parameters.
On 2016-Jan-21st ITRF2014 transformation parameter set was released. It has been added to Trafo_ITRF parameter set and the corresponding itrstrafo function has undergone some changes.
Just changed some minor things in the examples document
Just changed some minor things in the examples document
I was asked for a TM projection with latitude origin not being the equator (Luxembourgian LUREF).
Two bugs were found in helmert3D's parameter adjustment.
A mistyping was found in Helmert affine 3D transformation which is corerected now.
Small Bug in helmert3d due to the latest chenges "x0 is not found" .
Andrea Pinna of the University of Cagliari sent me improvements of the helmert- and helmertaffine-functions which allowed to save a lot of memory and time. I am glad he offered me to update my toolbox with his improvements. Thanks!
Fixed some minor warning issues.
Minor changes in the background document
Molodensky transformation added due to a user request. Minor bug in deg2dms fixed. In Matlab 2013b, helmert adjustments with worse normal equation matrices diverged while they converged in earlier versions. I changed that functions to overcome that.
There were two small bugs in helmertprojective3d.m which threw errors.
In itrstrafo.m there was a mistyping of two variable names which have been corrected. Additionally, a set of transformation parameters from IGS2005 to IGS2008 has been added in Trafo_ITRF.mat
A missing default was fixed in cart2ell.m
Update the coding error in d2affinetrafo.m Michael Schmitt mentioned on June 18th, 2012. Thanks, Michael!
Some minor changes in the existing functions were made. Additional transformation types have been added.
The existing files of the Geodetic toolbox hove been completely revised and many additional functionality has been added.
A left-out transpose-command rose an error in trafo3d with only two input parameters. Thanks to Jon Preston for noticing me.
There was some confusion with the additional file helmert3d which did not transform right-handed coordinate systems but left-handed, so the parameter sets did not fit in trafo3d.
I was asked to add a function which allows to calculate the transformation parameter set for cartesian datum transformation from two coordinate lists in both systems.
I was notified of a weak coded while condition which might cause the loop not to terminate although the criteria is already fulfilled. I changed that - results stay unchanged in any case.
I was asked by Matlab users without geodetic background to give additional information about the transformation steps and an example of usage. I therfore added a paper which explains the calculations and gives a short introduction.
The Projections.mat file now is also accessible for Matlab R13 users.
A big in adding the right zone identifier was found for usage with utm-coordinates when transforming from ellipsoidal to projection coordinates.
I forgot to enable backwards compability for the included mat-Files. Now the mat-Files also work with Matlab V6.x