|
interpolating mapping searching |
|
Inflection Point of a Great Circle:
The inflection point (apex) of a great circle is the point at which it achieves it's highest latitude. The antipode of this point is the point at which the great circle achieves it's lowest latitude. The latitude of these two points defines the latitude range of the great circle and one extreme of the latitude range of any great circle arc that contains one of them.
In Cartesian space we have three equations and three unknowns:
x2+ y2+ z2= r2 - the sphere.
ax + by + cz = 0 - the plane of the great circle.
z = zmax - the plane of the maximum latitude achieved by the great circle.
a, b, c and r are constants. x, y and z range over {-r, r}. The inflection point is where z is maximized so we need to find the (x, y) pair that maximizes z.
Solving for x and y we get:
x = (-cz - ((b2 * (-cz))/(a2 + b2)) (-+) (ab * sqrt(-(cz)2 +(r2 - z2)a2 +(r2 - z2)b2)/(a2+b2)))/a
and
y = (-bcz (+-) a * sqrt[-(cz)2 +(r2 - z2)a2 +(r2 - z2)b2]) / (a2 + b2)
Note the radical. If the term under the radical is negative there is no (real) solution. Consequently:
-(cz)2 +(r2 - z2)a2 +(r2 - z2)b2 >= 0
so
-c2z2 + r2a2 -z2a2 +r2b2 -z2b2 >= 0
r2a2 + r2b2 >= z2(c2+a2+b2)
(r2a2 + r2b2) / (c2+a2+b2) >= z2
which means
-sqrt[(r2a2 + r2b2) / (c2+a2+b2)]
<= z <=
+sqrt[(r2a2 + r2b2) / (c2+a2+b2)]
so
zmax = +sqrt[(r2a2 + r2b2) / (c2+a2+b2)]
which means
-(czmax)2 +(r2 - zmax2)a2 +(r2 - zmax2)b2
= -c2((r2a2 + r2b2) / (c2+a2+b2))
+ r2a2
-a2((r2a2 + r2b2) / (c2+a2+b2))
+ r2b2
-b2((r2a2 + r2b2) / (c2+a2+b2))
= -(a2 + b2 + c2)*((r2a2 + r2b2) / (c2+a2+b2)) + r2a2 + r2b2
= -(r2a2 + r2b2) + r2a2 + r2b2
= 0
and the inflection point is the point
((-c*zmax - ((-b2*c*zmax)/(a2+b2)))/a,
-(b*c*zmax)/(a2+b2),
zmax)
Determining whether a given point is inside or outside a given spherical polygon is rather difficult. Several effective algorithms exist for making the same determination on a plane, but those don't transfer well to the surface of the sphere. In the spheres package we have adapted the "point-at-a-distance" algorithm to work on the sphere for "reasonable" polygons, but it cannot be made absolutely foolproof.
The plane is infinitely large so finding a point outside a given polygon is relatively easy on the plane. One can simply pick (maximumX + 1, maximumY + 1) or, to be really safe, one can pick (maximumX * 10, maximumY * 10) where maximumX and maximumY are the maximum values of the given corner points. Indeed the algorithm is frequently refered to as the point-at-an-infinite-distance algorithm because developers just use the point (infinity, infinity), or substitute some absurdly large number for "infinity", no matter what polygon they are working with.
But the sphere is a closed surface so when you go far enough away from a given point you end up back where you started. Consequently there is no easy way to pick a point known to be external to the polygon. The spheres package employs a number of heuristics to guess an external point as outlined below, and in cases where that guess it likely to be wrong the SphericalPolygon class has a setExternalPoint() method to allow users to designate the known external point.
Given a point known to be external to the polygon determining if a given point is inside a given polygon is a simple matter.
Behavior is undefined if the given point is a corner point of, or on an edge of, the spherical polygon. In cases where the great circle arc constructed in step #1 just touches an edge of the polygon the algorithm may wrongly designate an external point as internal. But those cases should be excessively rare.
For most reasonable spherical polygons finding a known external point is fairly straightforward. A notable exception is a spherical polygon that is intended to outline an orbit of a satellite born sensor. The swath of a single orbit of a sensor on a polar orbiting satellite commonly covers both poles, all longitudes, and all latitudes while covering only a small portion of the Earth's surface. This is precisely why the spheres package includes an Orbit class, and for the purposes of working with spherical polygons it is assumed the polygon does not have these traits.
Another difference between polygons on the plane and polygons on the sphere is that while both merely separate the surface into two distinct areas there is no obvious way to define the "inside" of a spherical polygon. On the plane it is simple, the "inside" is the part with a finite area. But the sphere is a closed surface so both areas are finite. Moreover there is often good reason to define an area on the sphere that covers more than half the sphere.
Nevertheless the spheres package has to make some assumptions which are described below. Because the poles are discontinuities in the spherical coordinate system the spheres package generally avoids using them as an external point.
The spheres package uses the following heuristics to guess an external point.
This is clearly just a guess. One consequence of this guess is that if you are attempting to define an orbit as a polygon, and the swath of the sensor is wide enough to include both poles, the point-in-polygon algorithm will always be wrong because the "external" point is actually internal. If you are trying to use the point-in-polygon method with "unusual" polygons it is best to set your own external point. Polygons that include both poles count as "unusual".
Please send comments, questions and queries regarding the Spheres project to swick@geospatialmethods.org
Last modified: Mon Oct 14 18:24:51 2002