Geographical points + great circle = ellipse? (Matlab)

[MartinV]

I'm doing this in Matlab but it's not restricted to any particular software.

I have a bunch of geographical points (x,y coordinates for each) and I want to take all the points that are 50 km or closer to the reference point. I took the great-circle equation to convert geographical longitude and latitude into angular distance which I can then multiply with Earth's radius and I'm done.

The thing is when I plot all the points that are supposed to be 50 km away or closer to my center (and make sure both axis have the same scale by typing "axis equal"), the points make out an ellipse, not a circle. I rechecked the code and everything seems fine. I made double sure by drawing a circle on top of my plot and yes, some points stick out on the sides.

Do you guys have any idea what could cause this? I was thinking of maybe the data being distorted due to Earth's curvature but the data was collected from the Earth's surface, not from a satellite.


Here's my code that takes the events and sorts them according to their range:

function B = handBag(ref,A,R) %A is the main dataset, R radius, ref reference point

B = [];

for i = 1:length(A)
if (ang(ref,A(i, <= R/6371)
B = [B; A(i,:)];
end
end

end

function y = ang(A,B) %calculates the angle difference between two points

yy = sin(A(2)*pi/180) *sin(B(2)*pi/180) + cos(A(2)*pi/180) *cos(B(2)*pi/180) ...
*cos(-1*(A(1)-B(1))*pi/180);

y = acos(yy);
end

 

[Bill Simpson]

Create a small list of data points that you are ABSOLUTELY certain are all exactly 49 km from the center and are roughly uniformly positioned around the center. Use those with your software and see if they are all neatly a circle just inside your existing circle. Don't just "use your own code backwards" to generate these points because that could more easily just reproduce whatever errors you might already have, find some completely different independent way of getting these that you can be certain is correct.

If they are also an ellipse that tells you one thing. If some are inside and some are outside your circle that tells you something different.

Done right this should help you narrow down where the problem is by at least half.

 

[MartinV]

It turned out the problem was this:
https://en.wikipedia.org/wiki/Equirectangular_projection

The function I used to draw my circle on top of my points was this:

function [x,y] = circle(x0,y0,r)

alpha = 0:0.01:2*pi;
x = x0 + r/6371*180/pi *cos(alpha);
y = y0 + r/6371*180/pi *sin(alpha);
plot(x,y,'b-');
end

I changed it into this:

function [x,y] = circle(x0,y0,r)

alpha = 0:0.01:2*pi;
x = x0 + r/6371*180/pi *cos(alpha) /cosd(y0);
y = y0 + r/6371*180/pi *sin(alpha);
plot(x,y,'b-');
end