Geodetic to cartesian coordinates

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

Answer 1

[David Goodmanson]

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]

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]

Hi Benedict, After I realized why rotations are required I did quite a lot of revision to the answer above.

[Benedict Low]

Thanks David! Managed to get the pitch rotated.

[Tania]

Hi David,

I've been using the code you provided above to perform a similar rotational task.

In the code above you created a random array for team member data (11x2 arrary), assuming 11 rows of "team members" and one x,y coordinate for each team member.

I'm trying to rotate data in a 6135 x 15 table, where the data is layed out as follows 6135 data points of coordintes seperted into x & y columns for each person.

Time I Person 1 x I Person 1 y I Person 2 x I Person 2 y I Person 3 x I Person 3 y .....etc

When the data is graphed without rotation it looks like the below:

I would like to rotate this date to run length way along the Y axis:

Is there a way to adapt the code written above to rotate a large table of x,y coordinates, instead of single x,y data points?

Thnk you in advance for your help!

[David Goodmanson]

Hi Tanya,

yes there is. I'll assume that your table is converted to a 6135x15 matrix that I'll call At. Looks like there are 7 players. The 2x2 rotation matrix R works on an mx2 matrix of coordinates. The idea is to temporarily drop the time column to make 6135x14 matrix A , stack up vertically all the player x coordinates into a single 7*6135 column, do the same for y, concatenate those columns side-by-side, do the rotation with your R and then reconstruct a new 6135x14 matrix B. Then put the time back in.

A = At(:,2:end);              % temporarily drop the time
n = size(A,1)                 % times
p = size(A,2)/2               % players
x = A(:,1:2:end);             % n x p
y = A(:,2:2:end); 
x = x(:);                     % vertical stack, (n*p) x 1
y = y(:);         
xy = [x y];                   % (n*p) x 2
xyR = xy*R;                            
xR = reshape(xyR(:,1),n,p);   % n x p
yR = reshape(xyR(:,2),n,p);
B = zeros(size(A));           % n x (2*p)
B(:,1:2:end) = xR;
B(:,2:2:end) = yR;
Bt = [At(:,1) B];             % reinsert time column

Answer 2

[Meysam Mahooti]

%--------------------------------------------------------------------------

%

% Position: Position vector (r [m]) from geodetic coordinates

% (Longitude [rad], latitude [rad], altitude [m])

%

% Last modified: 2018/01/27 M. Mahooti

%

%--------------------------------------------------------------------------

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;

[Benedict Low]

Thank you Meysam! I will give it a try.