Prev: Matlab location function
Next: sign test, polarity
From: TideMan on 12 Mar 2010 19:01
On Mar 13, 11:34 am, "Roger Stafford"
<ellieandrogerxy...(a)mindspring.com.invalid> wrote:
> "Ha " <scifi...(a)126.com> wrote in message <hnd5rr$4e...(a)fred.mathworks.com>...
> > The story below happens on Earth surface.
>
> > There is an arc of a great circle (AB), the two endpoints of which are placed at (latA,lonA) and (latB, lonB). Beside, there is a quadrangle. The latitude and longitude of the four vertices (V1,V2,V3,V4) of the quadrangle are also known. The edges (e.g. V1V2,V2V4,...) of the quadrangle all lie on great cirles. Then, how can I know if the arc AB  crosses the quadrangle and calculate their intersects?
>
> -----------------
>   Your quadrangle question can be broken down into finding whether a pair of great circle arcs, each defined by a pair of points, will intersect, and if so what that point of intersection is.
>
>   I would suggest transforming your longitude and latitude coordinates (in radians) over to cartesian coordinates using the equations:
>
>  x = cos(theta)*cos(phi);
>  y = sin(theta)*cos(phi);
>  z = sin(phi);
>
> where theta is longitude measured positive going east of greenwich and negative west, and where phi is latitude measured positive going north from the equator and negative south.  We should have -pi<=theta<=pi and -pi/2<=phi<=pi/2.  These cartesian coordinates correspond to a hypothetical spherical "earth" of unit radius, but that does not interfere in the following computations.
>
>   The following is a sketch of how you might then proceed.  Let a, b, c, and d be vectors of the cartesian coordinate endpoints for the two arcs a-b and c-d obtained in this way.  Carry out the following computations:
>
>  p = cross(a,b); % p is along the normal to the plane of arc a-to-b
>  q = cross(c,d); % Similarly for q and arc c-to-d
>  t = cross(p,q); % t is along the line of intersection of the planes
>  s1 = dot(cross(p,a),t);
>  s2 = dot(cross(b,p),t);
>  s3 = dot(cross(q,c),t);
>  s4 = dot(cross(d,q),t);
>
>   You will find that if s1, s2, s3, and s4 are all of the same sign, the arcs will then and only then intersect.  In that case they intersect along +t if they are all positive and along -t if all are negative.
>
>   If they do intersect, you can transform the corresponding vector, t or -t, back into longitude and latitude (without worrying about its length.)  Letting x,y,z be its cartesian coordinates this reverse transformation can be accomplished this way:
>
>  theta = atan2(y,x);
>  phi = atan2(z,sqrt(x^2+y^2));
>
>   The above procedure ignores the indeterminate cases when 1) two ends of an arc are coincident or at antipodes (in which case p or q is of zero length) or 2) when the arcs lie in the same plane (in which case t is of zero length).  I will let you deal with these.
>
> Roger Stafford

Another way to do it is to use the gnomonic map projection to convert
the great circle lines to straight lines on a plane. Then the problem
becomes simple geometry.
You can do the transformation using the m_map toolbox (Google it).
 | 
Pages: 1
Prev: Matlab location function
Next: sign test, polarity