Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Geographical points + great circle = ellipse? (Matlab)

  1. Jul 26, 2013 #1
    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
     
  2. jcsd
  3. Jul 29, 2013 #2
    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.
     
  4. Jul 31, 2013 #3
    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
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook