# Geodetic to cartesian coordinates

34 views (last 30 days)
Benedict Low on 27 Mar 2017
Commented: Benedict Low on 1 Nov 2019
Hi,
I work with GPS units for footballers and I want to convert geodetic coordinates to cartesian coordinates. I have tried using Matlab's geodetic2enu function but the coordinates come out funny. As an example, here are the latitude and longitude values of each corner of a football pitch, converted using Matlab's " geodetic2enu" function.
Pitch_dimensions=[50.707025 4.206677;50.70665 4.205923;50.706062 4.207913;50.705676 4.207147];
[xEast, yNorth, zUp] = geodetic2enu(Pitch_dimensions(:,1),Pitch_dimensions(:,2),0,50.70665,4.205923, 0,wgs84Ellipsoid);
Pitch_dimensions=[xEast yNorth];
scatter(Pitch_dimensions(:, 1), Pitch_dimensions(:, 2));
The scatter plot should turn out to be a basic rectangle, but it shows a slanted pitch instead .
I just want to display a simple football field with a left-to-right direction of play, and the origin (0,0) to be the bottom-left corner flag of the pitch. Is there a rotation matrix that I need to use? Any help is appreciated.
Thank you.
Best regards,
Ben

David Goodmanson on 28 Mar 2017
Edited: David Goodmanson on 29 Mar 2017
Hello Benedict, (revised)
The four points do form a pretty good rectangle, but it doesn't look like it because the x and y axis scaling on the plot are different. Try 'axis equal' just after the plot command, which should show you a rectangle, hopefully. The following
pd = Pitch_dimensions % in meters
norm(pd(1,:)-pd(2,:)) % width
norm(pd(3,:)-pd(4,:)) % width
norm(pd(1,:)-pd(3,:)) % length
norm(pd(2,:)-pd(4,:)) % length
norm(pd(1,:)-pd(4,:)) % diag
norm(pd(2,:)-pd(3,:)) % diag
shows that there is a bit of inaccuracy in the coordinates, on the order of 1 m, but the diagonals agree within about 2 m so it's a rectangle. (I have to say that I don't have the geodetic2enu function so I had to make my own simplified version, but I think it works acceptably for this purpose).
Pretty long touch lines.
At this point I will assume that you have gone through the same geodetic2enu process for players as for the corners. The result is a 4x2 array for the corners and 11x2 arrays for the teams, plus maybe a 1x2 for the ball, all in meters. Some demo code for that is
% demo code for players and ball, unrotated, in meters
pd = Pitch_dimensions; % in meters
p42 = pd(4,:)-pd(2,:); % touch line
p12 = pd(1,:)-pd(2,:); % goal line
team1 = rand(11,1)*p42+ rand(11,1)*p12;
team2 = rand(11,1)*p42+ rand(11,1)*p12;
ball = rand*p42 + rand*p12;
Here is an example of rotation and plotting. I used plot instead of scatter because I think it is more versatile. Lots of possibilities for a plot. You will see that the boundary is a bit off due to 1 m inaccuracies but you could always make one manually.
% start of real code
pd = Pitch_dimensions;
bound = pd([2 4 3 1 2],:); % boundary points in ccw order for plotting
p42 = pd(4,:) - pd(2,:); % touch line
theta = atan2(p42(2),p42(1)); % theta is in radians
R = [cos(theta) -sin(theta); sin(theta) cos(theta)]; % rotation matrix
team1R = team1*R; team2R = team2*R; % rotation
ballR = ball*R; boundR = bound*R;
figure(1)
plot(boundR(:,1),boundR(:,2),'o-');
hold on
plot(team1R(:,1),team1R(:,2),'sb','MarkerFaceColor','b') % everton
plot(team2R(:,1),team2R(:,2),'sr','MarkerFaceColor','r') % liverpool
plot(ballR(:,1),ballR(:,2),'ok');
axis equal
hold off
I maybe went on too long but it's an interesting topic. Looking at the random player positions reminds me a lot of our intramural soccer team from college.

Benedict Low on 28 Mar 2017
Hi David,
Thank you for writing.
In addition to the pitch, I also need to convert the geodetic coordinates of 22 moving footballers to cartesian coordinates too. Thus, I am trying to find out, why the geodetic coordinates are slanted, and what can I do to adjust them.
Best regards, Ben
David Goodmanson on 29 Mar 2017
Hi Benedict, After I realized why rotations are required I did quite a lot of revision to the answer above.
Benedict Low on 31 Mar 2017
Thanks David! Managed to get the pitch rotated.

Meysam Mahooti on 1 Nov 2019
%--------------------------------------------------------------------------
%
% Position: Position vector (r [m]) from geodetic coordinates
%
%
%--------------------------------------------------------------------------
function r = Position(lon, lat, h)
R_equ = 6378.137e3; % Earth's radius [m]; WGS-84
f = 1/298.257223563; % Flattening; WGS-84
e2 = f*(2-f); % Square of eccentricity
CosLat = cos(lat); % (Co)sine of geodetic latitude
SinLat = sin(lat);
% Position vector
N = R_equ/sqrt(1-e2*SinLat*SinLat);
r(1) = (N+h)*CosLat*cos(lon);
r(2) = (N+h)*CosLat*sin(lon);
r(3) = ((1-e2)*N+h)*SinLat;

#### 1 Comment

Benedict Low on 1 Nov 2019
Thank you Meysam! I will give it a try.