New Mexico Tech
Earth and Environmental Science

ERTH 455 / GEOP 555 - Geodetic Methods

Lab 2: Coordinate Conversions

Note that you don't have to work on redoubt today!

Part I: Linear Algebra Review

Some useful Matlab functions:

• doc [function_name]: open MATLAB help (to function_name)
• ones(m,n): create mxn matrix of ones
• zeros(m.n): create mxn matrix of zeros
• eye(m,n): create mxn identity matrix
• diag(x): makes square diagonal matrix with the elements of x on diagonal.
• rand(m,n): create mxn matrix of random numbers
• dot(u,v): calculate dot product of vectors u, v
• rad2deg(x): convert radians to degrees
• deg2rad(x): convert degrees to radians

It's important to note that trig-functions in matlab expect their arguments to be in radians. E.g., sin(x) expects x to be in radians. Use sind(x) if x is in degrees (or, well, convert).

Vectors and Matrices:

• Making a row vector: x =[1 2 3 4 5]
• Making a column vector: x=[1;2;3;4;5] or x=[1 2 3 4 5]’
• The colon operator – incremental vectors, etc: x=[1:5] will give you [1 2 3 4 5] (increment by 1)
• Other increments: x=[0:0.25:2] gives [0 0.25 0.5 0.75 1 1.25 1.5 1.75 2] (this is useful in loops, too).
• Making a Matrix: A = [1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5] creates a 4x5 matrix
• Matrix from row vectors: x=[1 2 3 4 5]; A=[x;x;x;x] makes the same matrix as above
• Matrix from column vectors: x=[1 2 3 4 5]’; A=[x x x x] will give you a 5x4 matrix
• Matrix multiplication of matrices A, B: A*B
• Element-wise operations: A.*B, A./B, A.^2 operate on the individual elements of A (and B)

Writing a Function:

New functions can be saved to an .m and called in the same way as the built-in functions. For instance my_cool_function(a,b,c) would live in the file my_cool_function.m in the current directory or somewhere in the path.

To start out writing some_cool_function, you'll open the matlab editor, for instance:

    >>edit some_cool_function
    

which creates the new file in your current directory. Then you define the function. The generic way to do that is:

        function [r1,r2,...,rN] = some_cool_function(i1, i2, ..., iN) 
    

where r1,..,rN are individual output values, of which there can be none, and then the whole [...] = construct disappears, and i1,...,iN are input parameters that you'll use inside your function body.

Let's say you need a function that squares your input:

    >>edit square_it
    

and then you're adding this to the editor window:

    function squared = square_it(x)
        square = x.^2;
    end
    

Save it. Done. That's it. You can now use square_it on the matlab command line (if you're in the directory where square_it is saved):

    >> square_it(2)

    ans = 

          4

    >> square_it([2 3 4 5])

    ans = 

          4    9    16    25
    

That last line is the kicker: Note that I used the .^ operator to square x. This allows me to work on full vectors and square each element individually, rather than square a vector, which will not work! These vectorized operations are often quite efficient and you should try to make use of this whenever you can.

Exercise 1:

OK - let's get to work! How about some exercises to get your hands dirty on some linear algebra? Please turn in a script (*.m file) that I can run, which gives the results on the screen (i.e. omit ';' at the end of the lines)

• define 2 vectors u = [2;3]; v = [5 6]; and calculate the following:
• u*v
• v*u
• u'*v'
• Make them both column vectors and add them.
• The rotation matrix (in 2D) is defined as R = [cos(theta) -sin(theta); sin(theta) cos(theta)]
• note that for the following 2 exercises, it's likely a good idea to actually plot the vectors (starting at the origin) to check your results. The easiest way to do that is either using plotv(u) (requires the neural network toolbox to be installed) or plot([0, u(1)], [0, u(2)], ...).
• with u = [0.5, 0.5] rotate u by theta = 45, 90, 180 degrees by solving R*u (note that you'll have to recalculate R every time, or you could write R as a function that takes an angle and returns the rotation matrix; or write a loop that changes theta.
• If u=[0.5, 0.5] and v=[-0.6, 1.0] Using the dot-product and vector norms, what is the angle between u and v? Now rotate one of them far enough to make the dot product with the other vector zero. What does that mean?

Part II: Coordinate Conversion

We know that GPS coordinates are given in Earth-centered-Earth-fixed (ECEF) reference frame / coordinate system. However, those Cartesian coordinates are very hard to interpret and we need to be able to convert these back into ellipsoidal coordinates (and vice-versa). Today you'll be implementing 2 functions (in your favorite programming language) that do exactly this. You will need these in the following labs for your conversions.

For this exercise use the WGS84 ellipsoid definition:

semi-major a = 6378137.0 m
flattening f = 1/298.257223563
semi-minor b = -f*a+a

The handout from Hoffmann-Wellenhof, 2nd Ed., that I distributed as handout and PDF via canvas lays out the formulas. What you need to do is write some code that takes in arguments, applies the formulas, returns the results.

Exercise 2: Ellipsoid to ECEF coordinates

In your preferred language, write a function that takes WGS84 latitude, longitude, and height (above the ellipsoid!) and returns cartesian X,Y,Z coordinates. The formulas given in the handout are fairly straight forward. Note that you need to convert degrees to radians, though!

Exercise 3: ECEF coordinates to Ellipsoid

In your preferred language, write a function that takes ECEF cartesian coordinates X,Y,Z and returns WGS84 ellipsoid latitude, longitude, and height. I recommend to implement the iterative algorithm from the handouts! Again, note that you'll need to convert radians to degrees.

Test

Your code should convert between these ECEF and WGS84 coordinates for a GPS station in McMurdo, Antarctica station:

	ECEF		WGS84
X (m)	-1311703.172	lat (deg)	-77.8383494704
Y (m)	310814.982	lon (deg)	166.669329308
Z (m)	-6213255.16	height (m)	98.0221867058

The GPS station is MCMD. I grabbed the approximate XYZ coordinates from a receiver file (rinex) from 2014-12-31. This is by no means its precise position.

Deliverables: (submit via canvas!)

• A matlab script for part I, which prints all the requested results
• The 2 functions / scripts / programs that do the conversion from part II