Code covered by the BSD License  

Highlights from
Geodetic Toolbox

4.8

4.8 | 15 ratings Rate this file 169 Downloads (last 30 days) File Size: 42.7 KB File ID: #15285

Geodetic Toolbox

by

 

13 Jun 2007 (Updated )

Toolbox for angle, coordinate and date conversions and transformations. Version 2.97.

| Watch this File

File Information
Description

A collection of geodetic functions that solve a variety of problems in geodesy. Supports a wide range of common and user-defined reference ellipsoids. Most functions are vectorized. Most recent version can be found at <http://www.craymer.com/software/matlab/geodetic/>. Functions include:

Angle Conversions
  deg2rad - Degrees to radians
  dms2deg - Degrees,minutes,seconds to degrees
  dms2rad - Degrees,minutes,seconds to radians
  rad2deg - Radians to degrees
  rad2dms - Radians to degrees,minutes,seconds
  rad2sec - Radians to seconds
  sec2rad - Seconds to radians

Coordinate Conversions
  ell2utm - Ellipsoidal (lat,long) to UTM (N,E) coordinates
  ell2xyz - Ellipsoidal (lat,long) to Cartesian (x,y,z) coodinates
  sph2xyz - Shperical (az,va,dist) to Cartesian (x,y,z) coordinates
  xyz2sph - Cartesian (x,y,z) to spherical (az,va,dist) coordinates
  xyz2ell - Cartesian (x,y,z) to ellipsoidal (lat,long,ht) coordinates
  xyz2ell2 - xyz2ell with Bowring height formula
  xyz2ell3 - xyz2ell using complete Bowring version
  utm2ell - UTM (N,E) to ellipsoidal (lat,long) coordinates

Coordinate Transformations
  refell - Reference ellipsoid definition
  ellradii - Various radii of curvature
  ct2lg - Conventional terrestrial (ECEF) to local geodetic (NEU)
  dg2lg - Differences in Geodetic (lat,lon) to local geodetic (NEU)
  cct2clg - Conventional terrestrial to local geodetic cov. matrix
  clg2cct - Local geodetic to conventional terrestrial cov. matrix
  rotct2lg - Rotation matrix for conventional terrestrial to local geod.
  rotlg2ct - Rotation matrix for local geod. to conventional terrestrial
  lg2ct - Local geodetic (NEU) to conventional terrestrial (ECEF)
  lg2dg - Local geodetic (NEU) to differences in geodetic (lat,lon)
  direct - Direct geodetic problem (X1,Y1,Z1 + Az,VA,Dist to X2,Y2,Z2)
  inverse - Inverse geodetic problem (X1,Y1,Z1 + X2,Y2,Z2 to Az,VA,Dist)
  simil - Similarity transformation (translation,rotation,scale change)

Date Conversions
  cal2jd - Calendar date to Julian date
  dates - Converts between different date formats
  doy2jd - Year and day of year to Julian date
  gps2jd - GPS week & seconds of week to Julian date
  jd2cal - Julian date to calenar date
  jd2dow - Julian date to day of week
  jd2doy - Julian date to year & day of year
  jd2gps - Julian date to GPS week & seconds of week
  jd2mjd - Julian date to Modified Julian date
  jd2yr - Julian date to year & decimal year
  mjd2jd - Modified Julian date to Julian date
  yr2jd - Year & decimal year to Julian date

Error Ellipses
  errell2 - Computes error ellipse semi-axes and azimuth
  errell3 - Computes error ellipsoid semi-axes, azimuths, inclinations
  plterrel - Plots error ellipse for covariance matrix

Miscellaneous
  cart2euler- Converts Cartesian coordinate rotations to Euler pole rotation
  euler2cart- Converts Euler pole rotation to Cartesian coordinate rotations
  findfixed - Finds fixed station based on 3D covariance matrix
  pltnet - Plots network of points with labels

Example Scripts

  DirInv - Simple partial GUI script for direct and inverse problems
  DirProb - Example of direct problem
  Dist3D - Example to compute incremental 3D distances between points.
  InvProb - Example of inverse problem
  PltNetEl - Example plot of network error ellipses
  ToUTM - Example of conversion from latitude,longitude to UTM

Acknowledgements

This file inspired Utm2 Lat Long Geo Tiff Converter, Equinoxes And Solstices, and Geod Yet Another Geodetic Toolbox.

MATLAB release MATLAB 7.13 (R2011b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (41)
24 Sep 2014 Mike Craymer

giometar,
Unless otherwise stated in the help messages, you cannot use files of data as input to my toolbox functions. You must first read the data into Matlab variables and then use the variables as input to the functions.

In the example you give, if n.txt, e.txt and u.txt are text files containing local geodetic north, east and up coordinates, respectively, try using the following:

n=load('n.txt');
e=load('e.txt');
u=load('u.txt');
lat=0.76; % radians
lon=0.36; % radians
[dX,dY,dZ]=lg2ct(n,e,u,lat,lon);

Feel free to email me at mike@craymer.com if you need further help using my toolbox.

22 Sep 2014 giometar

>> [dx,dy,dz]=lg2ct('n.txt','e.txt','u.txt',0.76,0.36)
Error using .*
Matrix dimensions must agree.

Error in lg2ct (line 46)
dX=sum(RR'.*[dx dy dz],2);

28 Aug 2014 Mike Craymer

Christopher,
Please see my comment dated 27 Jul 2013 below where I reported the correction to the error in ell2utm. That comment also reports a correction to an error in utm2ell.

All,
Before reporting an error, please have a look at my comments for possible corrections. I was planning to post an update to the toolbox this summer with these and other fixes but, alas, have not been able to find the time. Hopefully in the next few weeks. I apologize for the delay.

27 Aug 2014 Christopher

In ell2utm.m, should
E2=lam.^3.*cos(lat).^3/6*(1-t.^2+h2);
be changed to:
E2=lam.^3.*cos(lat).^3/6.*(1-t.^2+h2);
?

I got an error as follows:
>> ell2utm(lon,lat)
Error using *
Inner matrix dimensions must agree.
Error in ell2utm (line 97)
E2=lam.^3.*cos(lat).^3/6*(1-t.^2+h2);

27 Aug 2014 Christopher

In ell2utm.m, should
E2=lam.^3.*cos(lat).^3/6*(1-t.^2+h2);
be changed to:
E2=lam.^3.*cos(lat).^3/6.*(1-t.^2+h2);
?

I got an error as follows:
>> ell2utm(lon,lat)
Error using *
Inner matrix dimensions must agree.
Error in ell2utm (line 97)
E2=lam.^3.*cos(lat).^3/6*(1-t.^2+h2);

03 Jul 2014 cai

大地主题正反算是关于测地线,而不是直线。
因此你的"direct"和"inverse"不是真正在解大地主题正反算。
"Geodetic problems" are about Geodesics on an ellipsoid, not straight line.
So your "direct" and "inverse" functions seems not working on "Geodetic problems".

12 Apr 2014 Mike Craymer

Many thanks to Kyle for reporting the error in utm2ell for points in the southern hemisphere. The function will be fixed in the next release following his suggestion to use negative zone values for such points. For now users should subtract 1e7 m from southern hemisphere northings.

10 Apr 2014 Kyle

Great work Mike.

Small problem with utm2ell I think though; it doesn't work for the southern hemisphere.

Specifically, this line does not modify the false northing for the southern hemisphere:

No(N>1e7)=1e7;

I would suggest perhaps using a signed Zone to denote southern hemisphere and then:

No(Zone<0)=1e7;

03 Mar 2014 Mike Craymer

Another error was discovered in the ell2utm function. The test for UTM zones less than zero should be a test for zones less than OR EQUAL TO zero. Replace

Zone=Zone+(Zone<0)*60-(Zone>60)*60;

with

Zone=Zone+(Zone<=0)*60-(Zone>60)*60;

The error occurs only for zone 60 (longitudes 174 - <180 east).

Thanks to Mohammad Ali Goudarzi for pointing out the error. I will post a new version of the toolbox with this correction and the ones given in my 27 Jul 2013 comment as soon as possible.

03 Mar 2014 Mike Craymer

Dear Dmitry,
Sorry I didn't see your post sooner. Here is an example of converting geocentric Cartesian coordinates to ellipsoidal with xyz2ell3 and then back again with ell2xyz:

>> [X,Y,Z]=ell2xyz(deg2rad(45),deg2rad(-79),300);
>> [lat3,lon3,h3]=xyz2ell3(X,Y,Z);
>> [X3,Y3,Z3]=ell2xyz(lat3,lon3,h3);
>> [X3-X;Y3-Y;Z3-Z]

ans =

1.0e-06 *

-0.1377
0.7087
0.7227

As you can see, the differences are less than a micron.

10 Sep 2013 Dmitry

Hi, thanks for the code!
I'm using the toolbox to compute movement across a custom size relatively flat ellipsoid (flattening = .8 for example). To get me a handle on the coordinate transformation functions, can someone please show an example of a surface point's cartesian coordinates being converted with xyz2ell3 and than back with ell2xyz. I cannot seem to arrive to a matching result. Thanks!

27 Aug 2013 Fei Liu  
27 Jul 2013 Mike Craymer

Two errors were discovered in the ell2utm.m and utm2ell.m functions. In ell2utm.m there is a vectorized calculation error for variable E2. Replace

E2=lam.^3.*cos(lat).^3/6*(1-t.^2+h2);

with

E2=lam.^3.*cos(lat).^3/6.*(1-t.^2+h2);

Thanks to John Marcovici for reporting this error.

And in utm2ell.m the function name on lines 1 & 2 is incorrect. Replace

function [lat,lon]=ell2utm(N,E,Zone,a,e2,lcm)
% ELL2UTM Converts UTM coordinates to ellipsoidal coordinates.

with

function [lat,lon]=utm2ell(N,E,Zone,a,e2,lcm)
% UTM2ELL Converts UTM coordinates to ellipsoidal coordinates.

Thanks to Andreas Wuestefeld for pointing out this error.

19 Apr 2013 wg

really amazing!!!

09 Jan 2013 Peng  
09 Aug 2012 Pavan

Helped a lot, Thanks.

25 Sep 2011 Mike Craymer

Re: 23 Jul 2011 post by muhammad faiz pa'suya

Dear Muhammad: The ell2utm function will automatically determine and output the standard UTM zone in which your input coordinates fall. To also have it output the central meridian for the zone(s), change the first line of the ell2utm.m function from

function [N,E,Zone]=ell2utm(lat,lon,a,e2,lcm)
to
function [N,E,Zone,lcm]=ell2utm(lat,lon,a,e2,lcm)

This will be implemented in the next version of the Toolbox. To understand how the zones are determined, see the Wikipedia entry for UTM.

22 Sep 2011 Shien Kwun Leong

Very useful! Thank you very much.

23 Jul 2011 joo tan

mike..
if i wnat to set the central meridian, like zon 50, where can i get the central meridian value??

19 Apr 2011 Vihang bhatt

It is really a good toolbox. Thank you for making my life little easy. One suggestion for the author on Julian dates. In the function jd2cal, jd has to be positive. jd can be negitive (if they exist), following modification needed to the code

if(a==0)
dy = c - e - fix(30.6001*f) + rem((jd+0.5),max(a,1));
else
dy = c - e - fix(30.6001*f) + rem((jd+0.5),a);
end

In above case you are able to get 1st gragorian date for zeroth julian day. otherwise greogrian calander start from 2nd jan.

Dr. Vihang Bhatt

13 Mar 2011 Mike Craymer

Sorry again! I keep forgetting that the MATLAB File Exchange doesn't like the use of the standard angle brackets enclosing URLs. Here is the clickable URL sans brackets:

http://www.craymer.com/software/matlab/geodetic/

13 Mar 2011 Mike Craymer

A new version 2.95 of the Geodetic Toolbox is now available. Unfortunately, I am unable upload it to the MATLAB Central File Exchange (keep getting an error). Until the Mathworks resolves the problem, you can find the new version at <http://www.craymer.com/software/matlab/geodetic/>.

Updates include:

- Added new jd2mjd & mjd2jd functions for modified Julian date conversions.
- Modified ct2lg & lg2ct functions to optionally use different lat,lon for each point.
- Corrected use of ell2utm function in ToUTM example script.
- Removed "clear all" in example scripts to avoid clearing user's workspace.
- Modified functions to use GRS80 reference ellipsoid by default.

See Content.m for the full list of changes.

21 Feb 2011 Mike Craymer

Sorry the URL's in my previous comment did not format correctly. They should have been specified as:

http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system

http://geographiclib.sourceforge.net/cgi-bin/GeoConvert

21 Feb 2011 Mike Craymer

Dear Ram,
By default, the ell2utm function uses the standard UTM zone definition where the earth is divided into 60 zones each 6 deg of longitude in width. Zone 1 is bounded by longtiudes 180 deg E and 186 deg E.

<http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system>

By this definition, Zone 44 is bounded by longitudes 78 deg E and 84 deg E, while Zone 45 is bounded by 84 deg E and 90 deg E. That is, Zone 44 starts at 180+(44-1)*6-360=78. The boundary between Zones 44 and 45 therefore bisects what you call Zone 45R. You can verify this with the on-line version of GeoConvert at

<http://geographiclib.sourceforge.net/cgi-bin/GeoConvert>

In any case, the ell2utm function allows you to specify the central meridian for any non-standard zone definitions. See help ell2utm for more information.

-Mike

21 Feb 2011 Mike Craymer

Dear Bryce,
You are right! The "clear all" in dates.m is not cool. I've replaced it with clearing the specific variables used in the script. A new version 2.95 of the toolbox will be posted shortly fixing this in dates.m as well as other affected example scripts. In the meantime, users should remove any "clear" or "clear all" when using the examples. Thank you for bringing this to my attention!
-Mike

07 Jan 2011 Ist

That's seems to be very useful package. Thank You!

30 Sep 2010 Ram Acharya

It is a good tool. However, I couldn't make its use while I wanted to convert lat longs to UTM, for lats [26.4546 -- 30.4485 degrees North ] and longs [80.074 -- 88.23319 degrees East]. The UTM zone is 45R . However your tool splits it into 44R and 45R. Could you review your code for that particular region. Thanks.

27 May 2010 Bryce

dates.m does a 'clear all' at the beginning of its loop. Not cool!

It would be much better to clear the actual names of variables that you use in the script. Resetting people's workspace variables is not nice.

Looks like a good toolbox, though. Thanks!

06 Feb 2010 Mike Craymer

Dear Jonathon,

I just saw your question about transforming lat,lon to N,E,U. The ct2lg function in your last step requires as input the X,Y,Z coordinate *differences* between the two points. In your example, you entered the X,Y,Z coordinates for just one point. ct2lg would have interpreted this as very large X,Y,Z coordinate differences. The correct conversion for your example is given below (you need to enter the X,Y,Z differences in ct2lg as Xec-Xec, Yec-Yec, Zec-Zec).

Thanks for your comments!
-Mike

>> [a,b,e2,finv] = refell('WGS84'); %get ellipsoid params
rlat = deg2rad(40); rlon = deg2rad(-108); %lat,lon
[Xec Yec Zec] = ell2xyz(rlat,rlon,0,a,e2) %convert to ECEF xyz
[N E Up] = ct2lg(Xec-Xec,Yec-Yec,Zec-Zec,rlat,rlon) %convert to local NEU

Xec =

-1.5119e+06

Yec =

-4.6532e+06

Zec =

4.0780e+06

N =

0

E =

0

Up =

0

16 Oct 2009 Jonathan

Hello,

Thank you for such a useful toolbox! I have a question, however, regarding transforming lat,lon on the WGS84 ellipsoid to local North East Up. If I use the same lat,lon coordinate as both my point to convert and my origin, why would I get a -21.5 km local North value? Shouldn't North and East both be zero?

Thanks again,
Jon

<SNIP>

>> [a,b,e2,finv] = refell('WGS84'); %get ellipsoid params
rlat = deg2rad(40); rlon = deg2rad(-108); %lat,lon
[Xec Yec Zec] = ell2xyz(rlat,rlon,0,a,e2) %convert to ECEF xyz
[N E Up] = ct2lg(Xec,Yec,Zec,rlat,rlon) %convert to local NEU

Xec =
-1.5119e+06
Yec =
-4.6532e+06
Zec =
4.0780e+06

N =
-2.1054e+04
E =
0
Up =
6.3693e+06

</SNIP>

29 Jul 2009 Sepideh

Thanks for that, very useful

24 Feb 2009 Mike Craymer

I corrected the bug for vector and array input to function ell2utm per Sebastian Holz's comments and posted a new version 2.91. Unfortunately, the fix only worked for row vectors. I've corrected it again in version 2.92 and tested it this time to make sure it works. Sorry for any problems this may have caused. In addition to error checking, I'll consider adding array input to all functions in the next major release (this Summer I hope).

18 Feb 2009 Sebastian Hölz

There is a potential bug in function ell2UTM, which may lead to wrong results in the calculated northing, if the funciton is called with a vector containing coordinates from both southern and northern hemisphere. Change lines 23-27, which read ...

if lat>=0
No=0;
else
No=10000000;
end

to ...

No(size(lat))=0; % False northing (north)
No(lat<0) = 1e7; % False northing (south)

Sebastian

27 Jan 2009 Christopher

under matlab which file do i put this in so that it is being used?

11 Sep 2008 jaime torres

gracias

09 Jul 2008 Gene Newton

Thank you.

01 Mar 2008 Ray Borough

Very Useful

25 Feb 2008 Jagkrich A.warangkul

Thanks

05 Oct 2007 Ali Ebrahimi  
28 Jun 2007 Nice Work

Nice work

15 Jun 2007 John D'Errico

This seems pretty good from what I checked. There are many conversions in here.

A couple of minor points. I notice that there is no error checking. For example, in ellradii one must provide an axis length and eccentricity (squared). But there are no checks on the proper number of parameters or even on the sign of the eccentricity.

In dms2deg, I noticed one interesting feature. Converting from 32 degrees, 12 minutes, -3 seconds, results in a degree prediction of:

dms2deg([32 12 -3])
ans =
-32.201

Overall, its still pretty good.

Updates
18 Feb 2009

Corrected ell2utm function for mixed north and south latitudes. Thanks to Sebastian Holz for pointing out the error and providing a fix.

24 Feb 2009

Corrected correction to ell2utm function in v2.91 (and tested it this time!).

31 May 2009

Added Matlab File Exchange BSD licensing. No changes to code since v2.92.

31 May 2009

Added Matlab File Exchange BSD licensing. No changes to code since v2.92.

11 Mar 2011

Added jd2mjd & mjd2jd for modified Julian date conversions. Modified ct2lg & lg2ct for different lat,lon for each point. Corrected ToUTM. Removed clear all in example scripts. Modified functions to use GRS80 reference ellipsoid by default.

31 Mar 2011
13 Feb 2013
13 Feb 2013

Added utm2ell, cart2euler & euler2cart; corrected cal2jd, ell2utm, xyz2ell; modified cct2clg & clg2cct; vectorized jd2yr; corrected help comments in lg2ct, pltnet. See individual functions for details.

Contact us