A collection of geodetic functions that solve a variety of problems in geodesy. Supports a wide range of common and userdefined 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 semiaxes and azimuth
errell3  Computes error ellipsoid semiaxes, 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
2.99.0.0  n/a 

2.99.0.0  Corrected function name & test for points in southern hemisphere in utm2ell; Corrected E2 computation, improved zone computation & corrected for zero Zones in ell2utm; Added transformations for both lefthanded LG & righthanded CT systems in simil. 

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

1.21.0.0  n/a 

1.18.0.0  n/a 

1.12.0.0  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. 

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

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

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

1.1.0.0  Corrected ell2utm function for mixed north and south latitudes. Thanks to Sebastian Holz for pointing out the error and providing a fix. 
Inspired: Antarctic geoid conversions, Equinoxes and Solstices, geod  yet another geodetic toolbox, UTM2LatLong GeoTiff Converter
Create scripts with code, output, and formatted text in a single executable document.
Ahmed Mohammed (view profile)
Mitchell Krouss (view profile)
Kana Nagai (view profile)
Rafay Ali (view profile)
Thank you very much. This is a very comprehensive list and very helpful.
Wei Dai (view profile)
ii's great ,very thanks
Wei Dai (view profile)
Hasan (view profile)
Great toolbox, exactly what I have been looking for! Thank you!
shufen Chen (view profile)
Hello Mike,I don't know how to use the function of geodetic toolbox to calculate the inclination of elliptic orbit, semimajor axis, eccentricity from satellite altitude, velocity, longitude and latitude. look forward to hearing from you!
Helen
Neil Gogoi (view profile)
Hamzeh Sadeghisorkhani (view profile)
Mike Craymer (view profile)
Dear Ankit,
I have already replied to you directly back in January but thought I should also post my reply to this forum for the benefit of others.
The r & v parameters in the dg2lg and lg2dg functions are the radii of curvature of the reference ellipsoid in the meridian and prime vertical, respectively. h is the ellipsoidal height (height above the reference ellipsoid) in the unit of length (usually meters).
The references I have used for the conversion between differences in geodetic latitude & longitude (dlat,dlon) and local geodetic (topocentric Cartesian) coordinates (dx,dy) are:
P. Vanicek & E.J. Krakiwsky. Geodesy: The Concepts. North Holland, 1986.
R.R. Sleeves. Mathematical Models for Use in the Readjustment of the North American Geodetic Networks. Technical Report No. 1, Geodetic Survey of Canada, Energy, Mines & Resources Canada, April 1984. [I can provide a PDF copy on request]
The meridian radius of curvature of the reference ellipsoid (v) is given by eqn. 1.4 in Steeves and eqn. 7.14 in Vanicek & Krakiwsky (denoted as M). The prime vertical radius of curvature of curvature (r) is given in the first line on page 20 of Steeves and by eqn. 15.58 in Vanicek & Krakiwsky (denoted as N). The formulae for converting local geodetic (dx, dy) coordinates to changes in latitude and longitude (dlat,dlon) are given in Steeves by eqns. 2.68 & 2.69. See also Fig 15.23 from Vanicek & Krakiwski illustrating the relationship between the Conventional Terrestrial (geocentric Cartesian) and local geodetic (topocentric Cartesian) coordinate systems.
Regards,
Mike
Andybai (view profile)
Ankit (view profile)
Hello Mike,
I have a question regarding the formulas that you have used in lg2dg and dg2lg functions. Can you provide me the reference for the formulas of v and r.
Also one more question, is the value for h in both the above mentioned functions equal to the longitude value (origin) in radians or something different?
Finally thanks a lot for sharing this toolbox. It helped me a lot. I am looking forward to hear back from you.
Regards,
Ankit
Mike Craymer (view profile)
Dear Denis,
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.
I have been meaning to upload a new version of the toolbox with the corrected utm2ell but, alas, never seem to find the time. Since my comment about the utm2ell error is now 3 years old and deeply buried amongst other comments, I have just posted the latest version 2.99 of my toolbox today. Sorry for taking so long to post it here.
Note that, as explained in the File Information notes, the latest version of the Geodetic Toolbox can always be found at http://www.craymer.com/software/matlab/geodetic/.
Denis Anikiev (view profile)
Thanks for the toolbox.
I've found that function utm2ell has wrong name inside (ell2utm) and also same typo in description
Mike Craymer (view profile)
Hi Amorim,
Sorry for the late reply. Been rather busy lately.
You can do what you want but you need to know the transformation parameters from the original system to the new one. Given the lat,lon of a point in the original system and an azimuth emanating from it, you can compute the azimuth in the new system as follows:
1) Solve the direct geodetic problem in the original system to compute the position of a second point (p2) an arbitrary distance away (say 1000 m) from the original point (p1). You can use my direct.m function for this.
2) Transform both the p1 & p2 points from the old system to the new system. If the transformation is a 7 parameter similarity (aka Helmert) transformation, you can use the latest version of my simil.m function that uses the geocentric Cartesian (CT) coordinate system. That version is not yet in the File Exchange but you can find it on my web site at
http://www.craymer.com/software/matlab/geodetic/
Note that you need to first convert the lat,lon of p1 & p2 to geocentric Cartesian coordinates to use simil.m. Use my ell2xyz.m for this conversion. The output transformed Cartesian coordinates of both points in the new system can then be converted to lat,lon using my xyz2ell.m.
3) With lat,lon of p1 & p2 in the new system, solve the inverse geodetic problem to get the azimuth and distance from p1 to p2. You can use my inverse.m function for this. The azimuth from inverse.m will be the azimuth in the new system you are looking for.
Rafhael Amorim (view profile)
Hi, Craymer.
This is a fantastic toolbox you are sharing with us.
I don't know, however, if it can help me with the following problem:
1) I have a geolocated point (latlong) and a path starting from this point to a particular given direction(bearing  measured from North Pole). How could I possibly translate this geodirection into a direction in other coordinate system?
Rob Whitfield (view profile)
Jaroslaw Tuszynski (view profile)
I needed some conversions and this toolbox was exactly what I needed. One confusing thing was NEU (NorthEastUp)coordinate system. I usually use NED (NorthEastDown) and I have seen other use ENU (EastNorthUp). But it is easy enough to flip the sign.
Mike Craymer (view profile)
Hi Jane,
That is correct. You can assign the output north (dg2lg's dx) and east (dg2lg's dy) to any variables to wish. Just make sure your y,x are not map projection coordinates. Topocentric Cartesian north and east are NOT the same as map projection northing and easting. If you want UTM map projection coordinates, use my ell2utm function.
You can, of course, define your local topocentric Cartesian system anyway you wish. However, the more common convention is to use x for north and y for east (lefthanded system), while map projections use y for north and x for east (righthanded system). Thus the warning about misinterpreting the output from dg2lg as map projection coordinates.
Feel free to contact me directly at mike@craymer.com if you want to discuss this further.
Jane (view profile)
Dear Mike,
Thank you for your reply and detailed explanations.
If I understood correctly, in a case where I want xaxis pointing East, yaxis pointing North, then I can still use your function dg2lg, but I need to assign the output as [y,x]. Am I right?
Thanks.
Jane
Mike Craymer (view profile)
Dear Jane,
The output of dg2lg is correct. The output dx,dy are in the local geodetic (LG) coordinate system centered at the lat,lon you specific in the input for dg2lg. The local geodetic coordinate system is defined by dg2lg, following the convention of Vanicek & Krakiwsky (Geodesy: The Concepts, 1986), as a lefthanded Cartesian coordinate system with the x axis pointing towards north (N), the y axis toward east (E) and the z axis up. This is defined in the help message for dg2lg.
In your example, dlon is larger than the dlat so that dy (E) will be larger than dx (N) as shown in the output from dg2lg. It appears the Woods Hole Institution's NDSF Utility is treating the dx as east and dy as north. That is the convention used for map projections like UTM, which the local geodetic system is not.
In a future version of my Geodetic Toolbox I will be changing the notation for the local geodetic system from dx,dy to dn,de to avoid the confusion with axes definition.
I hope this explains the dg2lg output for you. Don't hesitate to contact me if you have any further questions.
Jane (view profile)
Hi Mike,
Thanks a lot for sharing this file. I was using the function [dx,dy]=dg2lg(dlat,dlon,lat,lon). But the dx, dy results I got seems to be the dy, dx results I got from the website
http://www.whoi.edu/marine/ndsf/cgibin/NDSFutility.cgi?form=0&from=LatLon&to=XY
For example:
lat = 51.945086*pi/180;
lon = 4.039216*pi/180;
dlat = 51.947022*pi/180  lat;
dlon = 4.054263*pi/180  lon;
[dx,dy]=dg2lg(dlat,dlon,lat,lon)
dx =
215.4116
dy =
1.0347e+003
But the resutls from the above website are: dx = 1.0347e+003, dy = 215.4116.
May I know the possible reason for this?
Huarong
Val Schmidt (view profile)
OMG. Mike, you are totally correct. I was moving too quickly. How embarrassing. Thanks. Val
Mike Craymer (view profile)
Val,
Contrary to what you claim, the results from Palacios's deg2utm are indeed identical with my ell2utm function. You have apparently misinterpreted the output from deg2utm and ell2utm. In my ell2utm/utm2ell functions, n=Northing and e=Easting. These are explicitly defined in the help comments. In the notation of Palacios's deg2utm/utm2deg functions and many map projection formulae, x=Easting (what you call "a") and y=Northing (what you call "b"). Unfortunately, Palacios does not explicitly define what x & y are which may explain your misinterpretation. In the example you give, my ell2utm's Northing (n) and Easting (e) are identical with deg2utm's Northing (y, your b) and Easting (x, your a). Thus, ell2utm and deg2utm give identical results! I would therefore appreciate if you would reconsider your 3star rating of my toolbox. Thank you.
Val Schmidt (view profile)
I am comparing this code with deg2utm and utm2deg from Rafael Palacios available elsewhere in MATLAB central. I've been using these other functions for a long time, with other software and believe them to provide correct answers for conversion to WGS84 based UTM coordinates. However i don't get the same answers here and I'm not sure if I'm doing something wrong or it is a bug. The math in the implementations is sufficiently different that it's not easy for me to decipher. (This is not my science). Here's what I mean:
[a b c] = deg2utm(43,70)
a =
4.1849e+05
b =
4.7613e+06
c =
19 T
Compared with:
>> [a b e2 finv] = refell('WGS84')
a =
6378137
b =
6.3568e+06
e2 =
0.0067
finv =
298.2572
>> [n e z] = ell2utm(43*pi/180,70*pi/180,a,e2)
n =
4.7613e+06
e =
4.1849e+05
z =
19
This zone is correct, but the other values are vastly different. Any ideas what is amiss?
Thanks,
Val
Filip (view profile)
Very useful toolbox.
Geodetic network (1,2,3 D) adjustment toolbox would also be significant.
Filip (view profile)
Mike Craymer (view profile)
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.
giometar (view profile)
>> [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);
Mike Craymer (view profile)
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.
Christopher (view profile)
In ell2utm.m, should
E2=lam.^3.*cos(lat).^3/6*(1t.^2+h2);
be changed to:
E2=lam.^3.*cos(lat).^3/6.*(1t.^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*(1t.^2+h2);
Christopher (view profile)
In ell2utm.m, should
E2=lam.^3.*cos(lat).^3/6*(1t.^2+h2);
be changed to:
E2=lam.^3.*cos(lat).^3/6.*(1t.^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*(1t.^2+h2);
cai (view profile)
大地主题正反算是关于测地线，而不是直线。
因此你的"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".
Mike Craymer (view profile)
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.
Kyle (view profile)
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;
Mike Craymer (view profile)
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.
Mike Craymer (view profile)
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);
>> [X3X;Y3Y;Z3Z]
ans =
1.0e06 *
0.1377
0.7087
0.7227
As you can see, the differences are less than a micron.
Dmitry (view profile)
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!
Fei Liu (view profile)
Mike Craymer (view profile)
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*(1t.^2+h2);
with
E2=lam.^3.*cos(lat).^3/6.*(1t.^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.
wg (view profile)
really amazing！！！
Peng (view profile)
Pavan (view profile)
Helped a lot, Thanks.
Mike Craymer (view profile)
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.
Shien Kwun Leong (view profile)
Very useful! Thank you very much.
joo tan (view profile)
mike..
if i wnat to set the central meridian, like zon 50, where can i get the central meridian value??
Vihang bhatt (view profile)
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
Mike Craymer (view profile)
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/
Mike Craymer (view profile)
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.
Mike Craymer (view profile)
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/cgibin/GeoConvert
Mike Craymer (view profile)
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+(441)*6360=78. The boundary between Zones 44 and 45 therefore bisects what you call Zone 45R. You can verify this with the online version of GeoConvert at
<http://geographiclib.sourceforge.net/cgibin/GeoConvert>
In any case, the ell2utm function allows you to specify the central meridian for any nonstandard zone definitions. See help ell2utm for more information.
Mike
Mike Craymer (view profile)
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
Ist (view profile)
That's seems to be very useful package. Thank You!
Ram Acharya (view profile)
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.
Bryce (view profile)
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!
Mike Craymer (view profile)
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 XecXec, YecYec, ZecZec).
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(XecXec,YecYec,ZecZec,rlat,rlon) %convert to local NEU
Xec =
1.5119e+06
Yec =
4.6532e+06
Zec =
4.0780e+06
N =
0
E =
0
Up =
0
Jonathan (view profile)
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>
Sepideh (view profile)
Thanks for that, very useful
Mike Craymer (view profile)
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).
Sebastian Hölz (view profile)
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 2327, 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
Christopher (view profile)
under matlab which file do i put this in so that it is being used?
gracias
Thank you.
Very Useful
Thanks
Nice work
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.