Calculation of Cell Areas on an Oblate Planet

Purpose: To calculate the cell areas on an oblate spheroidal planet. A ``cell'' is defined as a region bounded by two meridians of longitude and two meridians of latitude. It is approximately a rectangular patch (except for cells near the two poles, which are spherical triangles). To better approximation, a cell is a trapezoid. Here, we rigorously derive the area.

We define f as the flattening of the sphere,

f =  a - b a

(1)

where a and b are the equatorial and polar radii, respectively. Note a > b for an oblate spheroid. For Earth, f = 0.00335364; for Mars, f = 0.00647630.

We begin with the definition of surface area in polar coordinates,

Area = ó õ l2 l1  ó õ b2 b1  r2 cosb  db dl

(2)

where l and b represent the longitude and latitude, respectively. For the sake of compactness, we do not carry the limits on the integrals in subsequent equations. However, the definite integrals remain to be evaluated after the dust settles.

The figure of an oblate spheroid is given (to order f²) by,

r = a æ è 1-fsin2b-  3 8 f2sin2 2b ö ø

(3)

Substituting (3) into (2) and doing the integration in the l direction gives:

Area = a2 Dl ó õ cosb æ è 1-fsin2b-  3 8 f2sin2 2b ö ø æ è 1-fsin2b-  3 8 f2sin2 2b ö ø   db

(4)

= a2 Dl ó õ cosb é ë 1 -fsin2b-  3 8 f2sin2 2b+f2sin4b+  3 8 f3sin2bsin2 2b -  3 8 f2sin2 2b+  3 8 f3sin2 2bsin2 2b+  9 64 f4sin4 2b ù û   db (5)

= a2 Dl ó õ é ë cosb -fsin2bcosb-  3 8 f2sin2 2bcosb+f2sin4bcosb+  3 8 f3sin2bsin2 2bcosb -  3 8 f2sin2 2bcosb+  3 8 f3sin2 2bsin2 2bcosb+  9 64 f4sin4 2bcosb ù û   db (6)

= a2 Dl ó õ é ë cosb- fsin2bcosb-  3 4 f2sin2 2bcosb+f2sin4bcosb+  3 4 f3sin2bsin2 2bcosb +  9 64 f4sin4 2bcosb ù û   db (7)

= a2 Dl é ë sinb -  f 3 sin3b-  3 4 f2 ó õ sin2 2bcosb  db+  f2 5 sin5b +  3 4 f3 ó õ sin2bsin2 2bcosb  db+  9 64 f4 ó õ sin4 2bcosb ù û   db (8)

= a2 Dl é ë I æ è sinb-  f 3 sin3b+  f2 5 sin5b ö ø   ê ê b2 b1  - II  3 4 f2 ó õ sin2 2bcosb  db   +  3 4 f3 ó õ sin2bsin2 2bcosb  db I   II+  9 64 f4 ó õ sin4 2bcosb  db IV   ù û (9)

To do the integrals in terms II-IV, note the following trigonometric identities:

sin2q = 2sinqcosq sin2 2q = 4cos2qcos4q = 4(sin2q- sin4q) (10a-10c) sin4 2q = 16(cos4q- 2cos6q+ cos8q) = 16(sin4q-2sin6q+ sin8q)

Term I is self-evident. Examining terms II-IV in turn,

II:

-  3 4 f2 ó õ æ è 16sin4 - 2bsin6b+ sin8b ö ø cosb  db = -  3 4 f2 é ë  16 5 sin5b-  2 7 sin7b+  1 9 sin9b ù û =-  12 5 f2sin5b+  3 14 f2sin7b-  1 12 f2sin9b ê ê b2 b1  (11)

III:

3 4 f3 ó õ sin2b×16 æ è sin4b- 2sin6b+sin8b ö ø cosb  db =12f3 ó õ æ è sin6b-2sin8b+sin10b ö ø cosb  db =12f3 é ë  1 7 sin7b-  2 9 sin9b+  1 11 sin11b ù û =  12 7 f3sin7b-  24 9 sin9b+  12 11 f3sin11b ê ê b2 b1  (12)

IV:

9 64 f4 ó õ 16 æ è sin4b-2sin6b+sin8b ö ø cosb  db =  9 4 f4 ó õ sin4bcosb  db-  9 2 f4 ó õ sin6bcosb  db+  9 64 f4 ó õ sin8bcosb  db =  9 4 f4  1 5 sin5b-  9 2 f4  1 7 sin7b+  9 64 f4  1 9 sin9b =  9 20 f4sin5b-  9 14 f4sin7b+  1 64 f4sin9b ê ê b2 b1  (13)

Putting it all together, Eq. (9) becomes:

Area = a2 Dl é ë I + II + III + IV ù û b2 b1

(14)

where b₁ and b₂ indicate the latitude limits to the region of integration.

Now time for the sanity checks:

·

First, in the case of a sphere, f=0, and all terms disappear except the lead term in I :

Areasphere = a2Dlsinb ê ê b2 b1

·

In this case, and for b₂=p/2 and b₁=-p/2 and Dl = 2p, we wind up with the area of a sphere, 4pa².

·

Also, in the limiting case of one of the poles, this expression reduces to a²DlDb/2, which is the expected area of the little triangle.

·

Second, in the limit of Dl® 0 and/or Db® 0, Area ® 0.

·

Third, all terms in the final expression are odd, indicating symmetry around the equator, as expected.

The following IDL procedure evaluates our expression for the case of Mars (a=3397 km):

`^^I= ^I

    `\`= `^^I=
  

SPHEROID_AREA


radius=3397.d0	; equatorial, in km


f=0.00647630	; flattening


gs=0.5d0	; gridsize, in degrees


halfgs=gs/2.d0


outfile='spheroid_area.dat'


	openw,unit,outfile,/get_lun


	printf,unit,'latitude Area(km2,sphere) Area(km2,spheroid)'


	for lat=90-halfgs,-90+halfgs,-gs do begin


	top=lat+halfgs	& bot=lat-halfgs


	area=!dtor*gs*radius2*(dotrig(top,f,gs)-dotrig(bot,f,gs))


	sphere_area=!dtor*gs*radius2*(sin(!dtor*top)-sin(!dtor*bot))


	printf,unit,lat,area,sphere_area,format='(1x,f6.2,2(3x,f12.6))'


endfor


	close,unit


	print,'Done. Data written to ',strn(outfile)


	DONE:


	stop


	end


function dotrig,lat,f,gs


	f2=f*f		&	f3=f2*f		&	f4=f3*f


	sinlat=sin(lat*!dtor)	&	sinlat2=sinlat*sinlat


	sinlat3=sinlat3	&	sinlat5=sinlat5


	sinlat7=sinlat7	&	sinlat9=sinlat9


	sinlat11=sinlat11


	term1=sinlat-(f/3.d0)*sinlat3+(f2/5.d0)*sinlat5


	term2=-(12.d0/5.d0)*f2*sinlat5+(13.d0/14.d0)*f2*sinlat7-(f2/12.d0)*sinlat9


	term3=(12.d0/7.d0)*f3*sinlat7-(8.d0/3.d0)*f3*sinlat9+(12.d0/11.d0)*f3*sinlat11


	term4=(9.d0/20.d0)*f4*sinlat5-(9.d0/14.d0)*f4*sinlat7+(f4/64.d0)*sinlat9


	print,term1,term2,term3,term4


	return,(term1+term2+term3+term4)


	end