Thread Subject:

how to judge if an arc of a great circle crosses a quadrangle and calculate the positions of the intersects of on a sphere surface

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?

Ha wrote:
> 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?

You would need some extremely detailed maps to determine that. Remember
that the Earth's surface varies quite irregularly near mountains and
near Death Valley.

The way your question was written hints that you are approximating the
Earth's surface as a sphere. That's not a very good approximation. You
can't even reasonably approximate "sea level" as a sphere -- at the very
least you need to use an oblate spheroid.

But calculations of the type you want to do are certainly much easier on
a sphere. Perhaps your instructor would be willing to modify the wording
of the problem?

Walter Roberson <roberson@hushmail.com> wrote in message <hnd6or$jlu$1@canopus.cc.umanitoba.ca>...
> Ha wrote:
> > 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?
>
> You would need some extremely detailed maps to determine that. Remember
> that the Earth's surface varies quite irregularly near mountains and
> near Death Valley.

Altitudes are also known. Oh, perhaps, I must take altitudes into consideration.

Does this work if A and B are not too far away?
d1 = distance([latA,lonA], [latB,lonB]);
d2 = sqrt(d1^2+(altA-altB)^2);


>
> The way your question was written hints that you are approximating the
> Earth's surface as a sphere. That's not a very good approximation. You
> can't even reasonably approximate "sea level" as a sphere -- at the very
> least you need to use an oblate spheroid.

Yes, you are correct that I consider it as a sphere. But, since the positions are relatively local, namely, the distances among the points are merely several degrees (up to 5 degrees), I thinks the discrepancy of distances on a sphere and a ellipse is trivil.

>
> But calculations of the type you want to do are certainly much easier on
> a sphere. Perhaps your instructor would be willing to modify the wording
> of the problem?

I have no instructor..... = =!

"Ha " <scifiles@126.com> wrote in message <hnd5rr$4eq$1@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

On Mar 13, 11:34 am, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> "Ha " <scifi...@126.com> wrote in message <hnd5rr$4e...@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).

"Ha " <scifiles@126.com> wrote in message <hnd5rr$4eq$1@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?
=================

Another approach. It might boil down to Roger's, but I'm not sure.
It assumes that the arc AB is not necessarily the shortest of the 2 arcs on the great circle between A and B, which you did not stipulate.

The points A and B are coplanar with the center of the earth (the coordinate origin) in a plane I'll denote pAB. For arc AB to intersect the quadrangle, it is equivalent to saying that plane pAB separates the 4 vertices of the quadrangle into separate half spaces. In other words if n is the normal to this plane, then

dot(n,V1), dot (n,V2), dot(n,V3), dot (n,V4),

are not all of the same sign. This provides a test of intersection.

Also, if two adjacent quadrangle vertices lie in different half spaces, you know an intersection occurs on the arc between those two vertices. Once you know between which two vertices an intersection occurs, solve for this intersection by solving for the intersection of

(1) plane, pAB
(2) the plane containing the two vertices and the origin, and
(3)the surface of the earth.

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hnem4s$b39$1@fred.mathworks.com>...
> Another approach. It might boil down to Roger's, but I'm not sure.
> It assumes that the arc AB is not necessarily the shortest of the 2 arcs on the great circle between A and B, which you did not stipulate.
> ........

  I made the fundamental assumption that all great circle arcs, arc AB and the four arcs of the quadrangle, are each the shortest distance on the surface of the sphere between their respective endpoints and therefore must be great circle arcs of no more than pi radians (180 degrees). That's crucial to the algorithm I presented. Otherwise it doesn't work.

  Of course if arcs greater than pi are allowed and if it is specified which of two possible arcs are chosen connecting each pair of vertices, then games could be played with the signs of the s quantities to find all possible intersection points. Without this additional information, however, the problem would not be well-defined. Even with this information it would be possible for two arcs to intersect in two antipodal points, leaving ambiguity as to which one ought to be chosen.

Roger Stafford

hi, guys. thanks for all your suggestions and discussions. I have looked throug and I'll read them carefully soon. All the points I'm considering are places in an area on the Earth surface with the longest distance up to ~5 degrees (~560km) and mean distance ~20km. Thus, only the smaller arc of a great circle is to be considered.

I have one other question: how can I take altitudes of the points into consideration? The difference of altitudes can be as large as 4km or even more, so altitude may be important.

"Ha " <scifiles@126.com> wrote in message <hnim5l$6ga$1@fred.mathworks.com>...

> I have one other question: how can I take altitudes of the points into consideration? The difference of altitudes can be as large as 4km or even more, so altitude may be important.
==============



In that case, what kind of curve would you want to use to connect the quadrangle vertices and the arc AB? With differing altitudes, it's obviously no longer a circular arc at the surface of a fixed-radius sphere. You could, for example, connect them via a straight line, or something else.

Also, what is the definition of intersection that you will use? If we're now connecting points with curves not confined to the surface of a sphere, these curves presumably have to intersect as curves in 3D space. What are the odds of that ever happening?

"Ha " <scifiles@126.com> wrote in message <hnim5l$6ga$1@fred.mathworks.com>...
> hi, guys. thanks for all your suggestions and discussions. I have looked throug and I'll read them carefully soon. All the points I'm considering are places in an area on the Earth surface with the longest distance up to ~5 degrees (~560km) and mean distance ~20km. Thus, only the smaller arc of a great circle is to be considered.
>
> I have one other question: how can I take altitudes of the points into consideration? The difference of altitudes can be as large as 4km or even more, so altitude may be important.
--------------
  Ha, I don't see that elevation or the elliptical shape of the earth enters into the problem as you have posed it. As I see it, your problem can be regarded as dealing only with angles, not the lengths of paths along the actual surface of the earth's irregular terrain. You can consider that all points have been radially projected with respect to the center of the earth onto a hypothetical perfectly spherical "earth", and you are, in effect, asking whether certain great circle arcs on this perfect sphere intersect and if so where that intersection is. You don't care what such an intersection's actual elevation is if you only ask for its longitude and latitude.

  Wikipedia, for example, defines longitude and latitude entirely in terms of angles and not elevations: "Latitude ... is the angle from a point on the Earth's surface to the equatorial plane, measured from the center of the sphere." "Longitude ... is the angle east or west of a reference meridian between the two geographical poles to another meridian that passes through an arbitrary point." "The combination of these two components specifies the position of any location on the planet, but does not consider altitude nor depth."

Roger Stafford

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hnin7d$mv0$1@fred.mathworks.com>...
> In that case, what kind of curve would you want to use to connect the quadrangle vertices and the arc AB? With differing altitudes, it's obviously no longer a circular arc at the surface of a fixed-radius sphere. You could, for example, connect them via a straight line, or something else.
>
> Also, what is the definition of intersection that you will use? If we're now connecting points with curves not confined to the surface of a sphere, these curves presumably have to intersect as curves in 3D space. What are the odds of that ever happening?


Oh, it's some confusing. I think I must state my question again.

At first, I deal with my problem on a spherical earth surface. I wanna find the intersection of an arc with a quandragle (in fact, I split a big quandragle into many small congruent grids and I have to solve many such intersection problems). All the lines involved are small arcs of great circle on the sphere. Then, I wanna calculate the length of each segment seperated by the small quandrangles, which is easily dealed with by using function distance(). However, Walter Roberson inspired me that altitude may be nontrival in the computation of distance for my problem. So, for simplicity, I still solve the intersection on a sphere, and then take the altitudes of all the points into consideration in the case of distance calculation. The altitude can be known by GPS records according to the lat and long. If two points have the same altitude, their distance is still calculated as a great
circle(thus the straight-line suggestion is invalid). But if one point is higher, their distance should be longer. I ever asked in the 3rd post if this works(whether it can be treated as a right triangle on a plane): d1 = distance([latA,lonA], [latB,lonB]); d2 = sqrt(d1^2+(altA-altB)^2);.

To conclude, my first problem is still the same, namely, the intersection-point problem on a sphere. And then, I wanna calculate the length of each segment with altitude involved. The maximum distance doesnot exceed 5 degrees.

(BTW, I'm dealing with the inversion of surface wave, actually.)

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hnj47l$7st$1@fred.mathworks.com>...
> Ha, I don't see that elevation or the elliptical shape of the earth enters into the problem as you have posed it. As I see it, your problem can be regarded as dealing only with angles, not the lengths of paths along the actual surface of the earth's irregular terrain. You can consider that all points have been radially projected with respect to the center of the earth onto a hypothetical perfectly spherical "earth", and you are, in effect, asking whether certain great circle arcs on this perfect sphere intersect and if so where that intersection is. You don't care what such an intersection's actual elevation is if you only ask for its longitude and latitude.

Hi, Roger. Yes, I don't care elevation on solving the intersection. See the post above for more details.

>
> Wikipedia, for example, defines longitude and latitude entirely in terms of angles and not elevations: "Latitude ... is the angle from a point on the Earth's surface to the equatorial plane, measured from the center of the sphere." "Longitude ... is the angle east or west of a reference meridian between the two geographical poles to another meridian that passes through an arbitrary point." "The combination of these two components specifies the position of any location on the planet, but does not consider altitude nor depth."
>
> Roger Stafford

Some complement:

For a small quandrangle V1V2V3V4, points V1 and V2 have the same latitude, so do points V3 and V4. Moreover, points V1 and V3 have the same longitude, so do points V2 and V4.

Points A and B are arbitrarily placed.

Ha wrote:

> Oh, it's some confusing. I think I must state my question again.
> At first, I deal with my problem on a spherical earth surface. I wanna
> find the intersection of an arc with a quandragle (in fact, I split a
> big quandragle into many small congruent grids and I have to solve many
> such intersection problems). All the lines involved are small arcs of
> great circle on the sphere.

Okay, that's an important clarification.

> Then, I wanna calculate the length of each
> segment seperated by the small quandrangles, which is easily dealed with
> by using function distance(). However, Walter Roberson inspired me that
> altitude may be nontrival in the computation of distance for my problem.

It was an accident ;-) What I was really trying to point out was that
when you take altitude into effect in joining the four points, then the
resulting curvalinear plane might not intersect a surface-level great
circle. To take a simple example, as Denver is "The Mile High City", the
planar surface of points nearby Denver would be about a mile above a
sea-level great circle arc. But putting everything on a uniform distance
sphere for the purposes of determining intersection will get around that
issue.

> But if one point is higher, their distance
> should be longer. I ever asked in the 3rd post if this works(whether it
> can be treated as a right triangle on a plane):

It isn't just whether one point is higher: if you are going to take
altitude into account, then the curved quad at altitude "sea level + X"
is bigger than the same angular coordinates at "sea level": the points
get further apart as the lat/long coordinates are projections towards
the center of the earth. Remember the classic formula for the length of
an arc: Pi*r*theta . For the purposes of determining intersection you
will use a fixed r, but for the purposes of real-world distances, you
will be using r that vary from point to point.

And no, it cannot be treated as a right triangle on a plane: the base
itself is curved and the upright is tangent to the curve of the earth at
that point, which is different than the tangent at the first point, so
your angles are not going to add up to 180.

You need to decide what _kind_ of distance you wish to measure. If you
use Euclidean distance, then for any two points on the surface of the
Earth, the shortest Euclidean distance is to go _through_ the Earth. If
you take the curvature of the Earth into account, then you have to
figure out what kind of arc you should use to connect two points of
different altitude. Tossing some ideas around in my head, I don't think
it would be unfair to use an arc length of Pi * theta * (r1+r2)/2, but I
will need to mentally double-check that idea.

 >Ha wrote
>> Oh, it's some confusing. I think I must state my question again.
>> At first, I deal with my problem on a spherical earth surface. I wanna
>> find the intersection of an arc with a quandragle (in fact, I split a
>> big quandragle into many small congruent grids and I have to solve
>> many such intersection problems). All the lines involved are small
>> arcs of great circle on the sphere.

You might also find some useful (non-altitude) hints at

http://www.darkshire.net/jhkim/programming/zipcodes/math.html

Walter Roberson <roberson@hushmail.com> wrote in message
> You might also find some useful (non-altitude) hints at
>
> http://www.darkshire.net/jhkim/programming/zipcodes/math.html

Thanks Walter. I've read it.

I have had some experiments and found that ellipticity and elevation are nontrival, which give rise to an error up to ~0.2%. An error less than 2% are absolutely acceptable for me since uncertainties in other measured data could be higher. Now I merely need to code the intersection problem. So, that's it.


Thanks to ALL :-) ! I'll return if I still have problem in the next work.