spherical geometry!
February 4, 2009 3:52 PM Subscribe
can you help me compute the area of a triangle on the sphere?
i have this pretty much down:
the formula for the area is:
a = R^2 ((A+B+C) - pi)
where:
a is area
R is the radius of the sphere
A, B and C are the triangle's angles
where i'm stuck is finding the angles A, B and C. i'm starting from 3 geographic (lat,lon) points. i think i pretty much get the principle there (finding the normals to the planes represented by the lines between the three points, etc.). but i can't quite make the leap.
wolfram and the rest have helped with the theory somewhat, but this math retard needs a bit of step-by-step hand-holding. i've worn out google trying to find a site that doesn't make a lot of unreasonable assumptions about my competence as a mathmatician.
i have this pretty much down:
the formula for the area is:
a = R^2 ((A+B+C) - pi)
where:
a is area
R is the radius of the sphere
A, B and C are the triangle's angles
where i'm stuck is finding the angles A, B and C. i'm starting from 3 geographic (lat,lon) points. i think i pretty much get the principle there (finding the normals to the planes represented by the lines between the three points, etc.). but i can't quite make the leap.
wolfram and the rest have helped with the theory somewhat, but this math retard needs a bit of step-by-step hand-holding. i've worn out google trying to find a site that doesn't make a lot of unreasonable assumptions about my competence as a mathmatician.
I'm no mathematician but, aren't A,B and C as angles always going to add up to 180? Shouldn't those be lengths as opposed to angles?
posted by Cat Pie Hurts at 4:27 PM on February 4, 2009
posted by Cat Pie Hurts at 4:27 PM on February 4, 2009
You've got latitude and longitude, so your spherical triangles are going to have to involve the North Pole somehow. Here's how I would tackle the problem:
Denote your points A, B and C in order of increasing longitude. Let N denote the North Pole. I'm going to assume for the sake of argument here that point B lies "inside" the triangle NAC (i.e., north of the great circle connecting A and C); see below for a similar argument if it lies outside. The area of the triangle ABC is then given by
Area(ABC) = Area(NAC) - Area(NAB) - Area (NBC)
For each of these three triangles on the right-hand side, you know the angle at the North Pole (given by differences in the longitudes) and the lengths of the two sides that meet at the North Pole (given by the latitudes of the two other points.) You can then use the law of cosines and the law of sines (as given on MathWorld) to figure out what the other two interior angles of each triangle are; and then you can use the area formula you mentioned above to figure out the area of each triangle.
If B lies south of the great circle connecting A and C, then you use roughly the same technique; the only difference is that the areas are related by
Area(ABC) = Area(NAB) + Area(NBC) - Area(NAC)
and you use the same procedure as above to calculate the area of each triangle.
Hope this helps. The process I've described is kind of involved, and perhaps some math wizard will pop in with some nicer identities that simplify it; but this'll get you to the answer, at least.
posted by Johnny Assay at 4:34 PM on February 4, 2009
Denote your points A, B and C in order of increasing longitude. Let N denote the North Pole. I'm going to assume for the sake of argument here that point B lies "inside" the triangle NAC (i.e., north of the great circle connecting A and C); see below for a similar argument if it lies outside. The area of the triangle ABC is then given by
Area(ABC) = Area(NAC) - Area(NAB) - Area (NBC)
For each of these three triangles on the right-hand side, you know the angle at the North Pole (given by differences in the longitudes) and the lengths of the two sides that meet at the North Pole (given by the latitudes of the two other points.) You can then use the law of cosines and the law of sines (as given on MathWorld) to figure out what the other two interior angles of each triangle are; and then you can use the area formula you mentioned above to figure out the area of each triangle.
If B lies south of the great circle connecting A and C, then you use roughly the same technique; the only difference is that the areas are related by
Area(ABC) = Area(NAB) + Area(NBC) - Area(NAC)
and you use the same procedure as above to calculate the area of each triangle.
Hope this helps. The process I've described is kind of involved, and perhaps some math wizard will pop in with some nicer identities that simplify it; but this'll get you to the answer, at least.
posted by Johnny Assay at 4:34 PM on February 4, 2009
The wikipedia page is actually pretty good. I think you convert the surface coordinates to spherical coordinates relative to the center of the sphere and then use those to determine tangent lines at the points. The angle between the tangent lines for each vertex of the triangle is the angle you want.
posted by GuyZero at 4:36 PM on February 4, 2009 [1 favorite]
posted by GuyZero at 4:36 PM on February 4, 2009 [1 favorite]
I'm no mathematician but, aren't A,B and C as angles always going to add up to 180?
Nope. You can have 3 right angles in a spherical triangle.
posted by GuyZero at 4:37 PM on February 4, 2009
Nope. You can have 3 right angles in a spherical triangle.
posted by GuyZero at 4:37 PM on February 4, 2009
Oh, I should mention that my procedure doesn't actually involve calculating the angles at A, B, and C directly; but you can find them from the other angles you derive in the procedure if you really want them.
I'm no mathematician but, aren't A,B and C as angles always going to add up to 180? Shouldn't those be lengths as opposed to angles?
Not on a sphere. Consider the triangle where one point is at the North Pole, one point is at 0° N/S, 0° E/W, and one point is at 0° N/S, 90 W. The interior angles of this triangle add up to 270°.
posted by Johnny Assay at 4:38 PM on February 4, 2009
I'm no mathematician but, aren't A,B and C as angles always going to add up to 180? Shouldn't those be lengths as opposed to angles?
Not on a sphere. Consider the triangle where one point is at the North Pole, one point is at 0° N/S, 0° E/W, and one point is at 0° N/S, 90 W. The interior angles of this triangle add up to 270°.
posted by Johnny Assay at 4:38 PM on February 4, 2009
Thanks, GuyZero. There's a reason that I'm no mathematician :)
posted by Cat Pie Hurts at 4:50 PM on February 4, 2009
posted by Cat Pie Hurts at 4:50 PM on February 4, 2009
Response by poster: GuyZero, that's basically where i'm at (angles betwen tangents are presumably the same as those between the normals of the planes formed by the great circles, which i was seeking in my question). the wiki article isn't a how-to though. i'm really stuck moving from that theoretical discussion to a practical application. can you point me to an example of the actual computation?
Johnny Assay, that's really clever. i'll try it out.
posted by klanawa at 5:45 PM on February 4, 2009
Johnny Assay, that's really clever. i'll try it out.
posted by klanawa at 5:45 PM on February 4, 2009
Write the points A,B,C in cartesian coordinates. Take the pair-wise cross-products of those three vectors; these are the normal vectors N_AB, N_AC, N_BC. Use the formula
v dot w = |v| |w| cos theta,
where v, w are vectors in R^3, and theta is the angle between v and w. This will give you the angles between the normal vectors. One can easily convince oneself that those angles are the same as the angles between the tangent lines. Then we apply the formula you quoted, and we are done.
However, the above method is not satisfactory, because the formula you quoted comes out of nowhere. An alternative solution would be the use the Gauss-Bonnet theorem. Or in our case Stokes theorem would suffice. To apply Stokes theorem, parametrize the boundary of the triangle. This is done, for instance, by finding the equation of the great circle passing through the pairs of points, parametrizing that circle in the usual way, then restricting the domain of the time variable. Then, do the appropriate line integral (I'll leave it to you to figure out which vector field F appearing in Stokes theorem to choose so that the double integral turns out to be the area).
posted by metastability at 6:28 PM on February 4, 2009 [2 favorites]
v dot w = |v| |w| cos theta,
where v, w are vectors in R^3, and theta is the angle between v and w. This will give you the angles between the normal vectors. One can easily convince oneself that those angles are the same as the angles between the tangent lines. Then we apply the formula you quoted, and we are done.
However, the above method is not satisfactory, because the formula you quoted comes out of nowhere. An alternative solution would be the use the Gauss-Bonnet theorem. Or in our case Stokes theorem would suffice. To apply Stokes theorem, parametrize the boundary of the triangle. This is done, for instance, by finding the equation of the great circle passing through the pairs of points, parametrizing that circle in the usual way, then restricting the domain of the time variable. Then, do the appropriate line integral (I'll leave it to you to figure out which vector field F appearing in Stokes theorem to choose so that the double integral turns out to be the area).
posted by metastability at 6:28 PM on February 4, 2009 [2 favorites]
I've forgotten all my math. But. Isn't there a way to find the total area of the surface of the sphere, and then divide that by the percentage of the surface of the sphere that the triangle-ish thing comprises? The three-dimensional version of finding the length of the arc on a circle (angle y) by multiplying the circumference by y/360?
Or is that a different way of saying what metastability said?
posted by gjc at 7:02 PM on February 4, 2009
Or is that a different way of saying what metastability said?
posted by gjc at 7:02 PM on February 4, 2009
Ugh, I feel compelled to hack away at this.
Another reference seems to indicate that you take the arc-length of each side (in degrees or radians) and then calculate the "spherical excess" E = sum of arc lengths - 180 deg.
Surface area A = R^2 * E.
If we take a right-right triangle that covers 1/8 of the sphere's surface we get
E = 90+90+90-180 = 90 or pi/2 radians
A = pi/2 * r^2 = (4pi*r^2)/8 or 1/8th of the surface of a sphere.
So I can prove it works in one case :)
Ah, a decent reference! The spherical excess measures the solid angle covered by the triangle in steradians and given the solid angle and the radius of the sphere, area is simple to calculate (solid angle * R^2).
posted by GuyZero at 9:50 PM on February 4, 2009
Another reference seems to indicate that you take the arc-length of each side (in degrees or radians) and then calculate the "spherical excess" E = sum of arc lengths - 180 deg.
Surface area A = R^2 * E.
If we take a right-right triangle that covers 1/8 of the sphere's surface we get
E = 90+90+90-180 = 90 or pi/2 radians
A = pi/2 * r^2 = (4pi*r^2)/8 or 1/8th of the surface of a sphere.
So I can prove it works in one case :)
Ah, a decent reference! The spherical excess measures the solid angle covered by the triangle in steradians and given the solid angle and the radius of the sphere, area is simple to calculate (solid angle * R^2).
posted by GuyZero at 9:50 PM on February 4, 2009
This other Wolfram reference may be of use. Mmmm... Haversines.
posted by GuyZero at 9:54 PM on February 4, 2009
posted by GuyZero at 9:54 PM on February 4, 2009
Yeah, based somewhat off of what metastability suggests, convert the spherical coords (lat & long) to cartestan coords and take the cross products to find the angles between the points.
Apparently there is no easy way to do cross-products in spherical coordinates so converting them to cartesian coordinates is the easiest way to do it.
so you have points a, b, c in spherical coordinates (same r value as they're on the same sphere)
calculate a', b', c' in Cartesian coordinates, ignoring the radius so they're unit length
calculate a' x b', b' x c' and a' x c'
take the inverse sine of each cross product to get the angle
add three angles, subtract pi radians
multiple by r^2 of sphere
I hope that will work. it seems easier than calculating the tangents.
posted by GuyZero at 10:11 PM on February 4, 2009
Apparently there is no easy way to do cross-products in spherical coordinates so converting them to cartesian coordinates is the easiest way to do it.
so you have points a, b, c in spherical coordinates (same r value as they're on the same sphere)
calculate a', b', c' in Cartesian coordinates, ignoring the radius so they're unit length
calculate a' x b', b' x c' and a' x c'
take the inverse sine of each cross product to get the angle
add three angles, subtract pi radians
multiple by r^2 of sphere
I hope that will work. it seems easier than calculating the tangents.
posted by GuyZero at 10:11 PM on February 4, 2009
Also, assuming the largest value you're getting back from inverse sine is pi, the biggest spherical triangle you can make is limited to being 2*pi*r^2, ergo any spherical triangle can be bounded by a hemisphere.
posted by GuyZero at 10:18 PM on February 4, 2009
posted by GuyZero at 10:18 PM on February 4, 2009
Upon further thought, what you really want to do is figure out the solid angle subtended between three vectors pointing from the center of the Earth to each of the three points. There's a formula on the Wikipedia page for "Solid Angle" (below the paragraph beginning, "An efficient algorithm for calculating the solid angle...") that looks like it would work perfectly. Just convert the latitudes and longitudes to Cartesian coordinates, take the appropriate dot products and triple products, and plug and chug.
posted by Johnny Assay at 5:44 AM on February 5, 2009
posted by Johnny Assay at 5:44 AM on February 5, 2009
Also, having thought about this WAY too much, I have basically rediscovered the formula you gave in the first place but have switched interior angles for the sides for the actual corners of the triangle.
Second, a set of three points on a sphere actually describes more than one triangle. So the result you're going to get from this formula is, I believe, the area of the smallest triangle, though I'm not sure if I can prove that. The issue is that using the cross product to compute the angle always gives you the interior angle (I think) when for any given side you could have used the interior angle or the exterior angle. (This may only be true in pairs of angles, but whatever)
posted by GuyZero at 9:23 AM on February 5, 2009
Second, a set of three points on a sphere actually describes more than one triangle. So the result you're going to get from this formula is, I believe, the area of the smallest triangle, though I'm not sure if I can prove that. The issue is that using the cross product to compute the angle always gives you the interior angle (I think) when for any given side you could have used the interior angle or the exterior angle. (This may only be true in pairs of angles, but whatever)
posted by GuyZero at 9:23 AM on February 5, 2009
Define the vectors r1, r2, and r3 as the vectors pointing from the center of the sphere to each respective point, divided by the radius of the sphere R (so that all are unit vectors), and r12 as the vector tangent to the great circle connecting point 1 to point 2 at point 1 (note that the order of 1 and 2 is important in the definition). What we are looking for is the angle between r12 and r13, for example.
The great circle connecting points 1 and 2 is contained within a plane defined by the vectors r1 and r2. A normal to this plane is given by (r1xr2). The vector r12 should be perpendicular to both r1 and (r1xr2). Up to a possible sign, it is given by r1x(r1xr2).
The angle t23 between r12 and r13 is then given by:
Cos[t23] = (r1x(r1xr2)).(r1x(r1xr3))
A little simplification using the triple product rule and the definition of a unit vector gives:
Cos[t23] = (r2.r3) - (r1.r2)*(r1.r3)
The remaining angles are generated by permuting the labels 1, 2, and 3:
Cos[t13] = (r1.r3) - (r2.r1)*(r2.r3)
Cos[t12] = (r1.r2) - (r1.r3)*(r2.r3)
You can verify that the dot products are given by:
r1.r2 = Cos[p1 - p2] Cos[t1] Cos[t2] + Sin[t1] Sin[t2]
where t1, t2 are the latitude and p1, p2 are the longitude of points 1 and 2, respectively.
posted by dsword at 9:44 AM on February 5, 2009
The great circle connecting points 1 and 2 is contained within a plane defined by the vectors r1 and r2. A normal to this plane is given by (r1xr2). The vector r12 should be perpendicular to both r1 and (r1xr2). Up to a possible sign, it is given by r1x(r1xr2).
The angle t23 between r12 and r13 is then given by:
Cos[t23] = (r1x(r1xr2)).(r1x(r1xr3))
A little simplification using the triple product rule and the definition of a unit vector gives:
Cos[t23] = (r2.r3) - (r1.r2)*(r1.r3)
The remaining angles are generated by permuting the labels 1, 2, and 3:
Cos[t13] = (r1.r3) - (r2.r1)*(r2.r3)
Cos[t12] = (r1.r2) - (r1.r3)*(r2.r3)
You can verify that the dot products are given by:
r1.r2 = Cos[p1 - p2] Cos[t1] Cos[t2] + Sin[t1] Sin[t2]
where t1, t2 are the latitude and p1, p2 are the longitude of points 1 and 2, respectively.
posted by dsword at 9:44 AM on February 5, 2009
Response by poster: so, dsword, when you say this:
posted by klanawa at 8:20 PM on February 5, 2009
A normal to this plane is given by (r1xr2)you are taking the cross product of r1 and r2, right? so what is r1? as above, finding the cross products of spherical coordinates isn't easy, so are these cartesian points? if so, is this reasonable?
x = r1.y*r2.z - r1.z*r2.y;
y = r1.z*r2.x - r1.x*r2.z;
z = r1.x*r2.y - r1.y*r2.x;
i think this is my primary hangup, but as you can see, i'm still a little lost. i've integrated the advice into a solution, but the output seems to be wrong, so i'm trying to track down my error.posted by klanawa at 8:20 PM on February 5, 2009
Yes, I mean the cross product, but you don't need to actually calculate any cross products. The expressions simplify to ones involving only dot products because I am dotting triple products together, and we have the triple product rule:
Ax(BxC) = (A.C)B - (A.B)C
r1 points from the center of the sphere to the point 1 on the surface of the sphere. It is unit-normalized (i.e., in spherical coords, it is given by (1,t1,p1), in Cartesian coords, it is given by (Cos[t1]*Cos[p1], Cos[t1]*Sin[p1], Sin[t1]).
But again, you don't need any of this. To find the angle at vertex 1, just calculate each of the dot products that you need:
(r1.r2) = Cos[p1 - p2] Cos[t1] Cos[t2] + Sin[t1] Sin[t2]
(r1.r3) = Cos[p1 - p3] Cos[t1] Cos[t3] + Sin[t1] Sin[t3]
(r2.r3) = Cos[p2 - p3] Cos[t2] Cos[t3] + Sin[t2] Sin[t3]
(these have been calculated by converting to Cartesian coords and are given directly by the quantities you know--latitude (t1,t2,t3) and longitude (p1,p2,p3)), and plug them into the appropriate formula, in this case it is:
Cos[t23] = (r2.r3) - (r1.r2)*(r1.r3)
Take the arccosine of both sides to find t23.
Calculating the Cartesian components is an intermediate step, which I did to derive the expressions for (r1.r2), etc... There's no need for you to go through that.
posted by dsword at 10:18 AM on February 6, 2009
Ax(BxC) = (A.C)B - (A.B)C
r1 points from the center of the sphere to the point 1 on the surface of the sphere. It is unit-normalized (i.e., in spherical coords, it is given by (1,t1,p1), in Cartesian coords, it is given by (Cos[t1]*Cos[p1], Cos[t1]*Sin[p1], Sin[t1]).
But again, you don't need any of this. To find the angle at vertex 1, just calculate each of the dot products that you need:
(r1.r2) = Cos[p1 - p2] Cos[t1] Cos[t2] + Sin[t1] Sin[t2]
(r1.r3) = Cos[p1 - p3] Cos[t1] Cos[t3] + Sin[t1] Sin[t3]
(r2.r3) = Cos[p2 - p3] Cos[t2] Cos[t3] + Sin[t2] Sin[t3]
(these have been calculated by converting to Cartesian coords and are given directly by the quantities you know--latitude (t1,t2,t3) and longitude (p1,p2,p3)), and plug them into the appropriate formula, in this case it is:
Cos[t23] = (r2.r3) - (r1.r2)*(r1.r3)
Take the arccosine of both sides to find t23.
Calculating the Cartesian components is an intermediate step, which I did to derive the expressions for (r1.r2), etc... There's no need for you to go through that.
posted by dsword at 10:18 AM on February 6, 2009
OOPS! I forgot to renormalize the cross-product. It's just an overall factor difference:
Cos[t23] = ((r2.r3) - (r1.r2)*(r1.r3)) / ( Sqrt[1-(r1.r2)^2]*Sqrt[1-(r1.r3)^2)] )
Cos[t12] = ((r1.r2) - (r1.r3)*(r2.r3)) / ( Sqrt[1-(r1.r3)^2]*Sqrt[1-(r2.r3)^2)] )
Cos[t13] = ((r1.r3) - (r1.r2)*(r2.r3)) / ( Sqrt[1-(r1.r2)^2]*Sqrt[1-(r2.r3)^2)] )
That should work much better. In particular, you should get that the interior angles of small triangles add up to pi.
posted by dsword at 12:24 PM on February 6, 2009
Cos[t23] = ((r2.r3) - (r1.r2)*(r1.r3)) / ( Sqrt[1-(r1.r2)^2]*Sqrt[1-(r1.r3)^2)] )
Cos[t12] = ((r1.r2) - (r1.r3)*(r2.r3)) / ( Sqrt[1-(r1.r3)^2]*Sqrt[1-(r2.r3)^2)] )
Cos[t13] = ((r1.r3) - (r1.r2)*(r2.r3)) / ( Sqrt[1-(r1.r2)^2]*Sqrt[1-(r2.r3)^2)] )
That should work much better. In particular, you should get that the interior angles of small triangles add up to pi.
posted by dsword at 12:24 PM on February 6, 2009
« Older When your lawyer malpractices you (remember to... | Need help hunting down a parasitic worm contracted... Newer »
This thread is closed to new comments.
What I'd do to find angle A involves a great circle sailing to get the initial course A-B, then one to get A-C, then a bit of squinting and muttering to work out what to do with those numbers to get that angle; it'd probably be easier to draw a diagram. The point being that the numbers you get from this handy great circle calculator will all come out as compass courses, whereas if you do it by hand you'll get quadrantal courses, which can get a bit more confusing. Yu could probably just subtract the smaller one from the larger one, in most cases....
If the above, which I think sums up as 'here is someone else's computer doing most of the maths for you', is enough, good luck. If you'd like any more explanation with real maths in, hopefully someone else will be along shortly, but otherwise I might be more use when sober, I'll have another look then.
posted by Lebannen at 4:21 PM on February 4, 2009