"good 'nuff" approximation

[jim mcnamara]

We calculate bounded areas - the boundaries are given in latitude/longitude ie, 48.5 deg N to 48.8 deg N, 105 deg W to 105.2 deg W. The bounded area is approximately the area of a trapezoid, except at the North Pole. The problem is 'approximately', as it gets more/less approximate as latitude changes. If the quadrilateral is small enough and not too far North,we just pretend it's square and go on.

The problem comes for larger ones or one in the far North.

Breaking the large area up into really short trapezoids, gets more accurate results. No surprise. We need accuracy circa 1:10000.
I can test for that. Run several tests using smaller trapezoids until we get where we want to be, or results do not change. This means a lot of extra cpu cycles.

The problem is performance. We do several million of these a day, and they are part of real-time data acquistion. I cannot waste time narrowing down the best approximation. I need to jump in at the reasonable level from the get-go. Right now we have a bunch of data in tables and have lookups that return heuristics. Works. Too much overhead (i/o on remote Oracle kernel is too slow)

Here is some dummy sample data to show you:
The red lines are the point where we got reasonably close to a decent result. They good result point can be different, but you can't tell a priori where it will be. I can't anyway.

from 0 to 5, 0.10 wide
stepsize:1.00000 area:34389.4801
stepsize:0.10000 area:35075.3246
stepsize:0.01000 area:34459.6559
stepsize:0.00100 area:34391.2357
stepsize:0.00010 area:34391.2357

from 5 to 10, 0.10 wide
stepsize:1.00000 area:33856.3562
stepsize:0.10000 area:34526.3349
stepsize:0.01000 area:33924.9792
stepsize:0.00100 area:33864.8234
stepsize:0.00010 area:33858.8076

from 15 to 20, 0.10 wide
stepsize:1.00000 area:31307.5109
stepsize:0.10000 area:31308.4049
stepsize:0.01000 area:31308.4238
stepsize:0.00100 area:31308.4240
stepsize:0.00010 area:31309.0321

from 20 to 25, 0.10 wide
stepsize:1.00000 area:29367.1160
stepsize:0.10000 area:29369.1078
stepsize:0.01000 area:29369.1277
stepsize:0.00100 area:29369.1279
stepsize:0.00010 area:29369.6933

from 85 to 90, 0.10 wide
stepsize:1.00000 area:75.1336
stepsize:0.10000 area:78.6199
stepsize:0.01000 area:78.6441
stepsize:0.00100 area:78.6444
stepsize:0.00010 area:78.6444

As you can see, it's not necessarily predictable, as far as I can tell.
Does anyone know of a non-exhaustive approach that will provide
a quick, good guess as to the best choice for "stepsize" - the height of the small trapezoids?

Human-friendly versions of the two c modules that do this are embedded in a file i/o routine to make the above junk data. It's attached. The code actually stays memory-resident on a unix box, gets data from all over, talks to Oracle, generates areas, and sends it to somebody else. Did not include all that.

And the actual modules use inlining & ASM. Of course, I've profiled them for best performance.

[Guv]

I have forgotten the details and have none of my texts on the subject.

There is a discipline called Differential Geometry. That discipline describes a formula called The First Fundamental Form, which provides an function which can be integrated to yield the arc length along the surface. Three sub formulae from this Form can be used to specify a function which can be integerated to determine area on the surface.

Both of the above can be used to make approximating formulae for arc length and surface area.

I have not done a good job of describing the above, and it probably does not convey anything useful to you.

My purpose is to supply you with some jargon for use in a Web search. Better yet, find a college mathematics major and feed him the jargon. He will be able to provide you with an approximating formula for the areas you are working with. Just tell him you are not intereseted in theory and just want a formula to approximate the area.

[NotLKH]

So, what you need is a formula for calculating the surface area of a curved trapezoid defined by Lat/Long Values on a sphere?

I've forgotten my Calc a LONG time ago, but isn't this an answerable calc question?

[Guv]

Just found an old text of mine. The following is an approximating formula for area on a sphere.

Area = R^2*cos(Latitude)*DeltaLatitude*DeltaLongitude

Where Latitude is zero at equator and 90 degrees at the North Pole.

The above refers to a patch of area bounded by two parallels of Latitude and two Longitude lines. An approximately trapezoidal area.

R is the radius of the sphere (Earth, I suppose). Note that the Earth is not exactly a sphere. The Equatorial radius is approximately 3964 miles, and the Polar radius is approximately 3950.5 miles. If you use miles for radius, area will be in square miles. If kilometers, then square kilometers.

DeltaLatitude is the difference (radians) between the two parallels of Latitude.

DeltaLongitude is the difference (radians) in Longitude.

Using average latitude for cos function should be reasonable.

It seems reasonable to break a big area into a lot of little ones. You could manage to use the same DeltaLatitude and DeltaLongitude for many small areas.

If you divided a large trapezoid-like area into a 10X10 checker board of 100 smaller ones, you could use a formula like the following (ten not factored out to make the formula a bit more understandable).

R^2*[ 10*cos(Lat1) + 10*cos(Lat2) + . . . + 10*cos(Lat10) ]*DeltaLat*DeltaLong

For computational purposes, the ten would be factored out.

Perhaps you should find somebody to verify the above. I think it is reasonable. Pehaps you could do some experimental calculations to verify formula and determine precision obtained for various size trapezoidal areas.

[jim mcnamara]

We'll try Guv's version. 10 ^-4 the CPU usage.

Thank you.

[Guv]

I did not pay enough attention to the original post.

The desired height of the trapezoids should be a function of latitude only. It seems to me that you need smaller increments as you get near the North Pole.

You should be able to experiment and build a table which shows the best increment for various ranges of latitude.

The formula I posted from a differential geometry text is correct for small areas. It treats an area as a small rectangle. Modifying it to treat the area as a small trapezoid would improve it. Also, I made a mistake in thinking that a checker board approximation was a good idea. Breaking a long East-West strip into smaller areas does not improve accuracy.

Consider the following formula.

Area = Radius^2*cos(Latitude)*DeltaLatitude*DeltaLongitude.

DeltaLatitude is the height of an approximate rectangle in radians, while DeltaLongtitude is the length. It would probably improve accuracy to make a trapezoid formula (such a formula is probably implemented by your C Code).

Area = R^2*[ cos(Lat1) + cos(Lat2) ]*DeltaLat*DeltaLong / 2

The units for radius determine the units for the area.

The cosine value is the critical part of the formula. Following are cosines for various Latitudes.

00 ---- 1.000 000
01 ----- .999 848
02 ----- .999 391
03 ----- .998 630

44 ----- .719 340
45 ----- .707 107
46 ----- .694 658

87 ----- .052 336
88 ----- .034 899
89 ----- .017 452

The percent change is obviously greatest near the pole. There is where you need the thinnest trapezoids. A little bit of experimentation should provide the data you need to assign height increments to ranges of latitude.