# What projection to use to calculate polygon areas of city to global size from lat/lon data?

I have a set of bounding boxes in WGS:84, and I would like to calculate their areas using Python.

• The bounding boxes are of all sizes, from city sized to covering the entire Earth.
• The bounding boxes may cross the dateline, but I'd imagine >99% of them don't.
• The areas don't need to be exact. To within 100^2 km is fine.

I will be doing the reprojection using pyproj in Python, and the area calculation's probably using Shapely. The question is - what projection should I use?

Wikipedia has an entire set of equal-area projections - https://en.wikipedia.org/wiki/Map_projection#Equal-area

Over on SO, there's a question dealing with this, but the answers all seem to have caveats that make them unsuitable for my purpose:

Which projection do I reproject to before doing the area calculation?

Edit

If I use the maths from here which @barrycarter helpfully provided, I end up with this Python:

``````import math

# Radius of earth in km^2
R = 6371
lat1 = -90
lat2 = 90
lon1 = -180
lon2 = 180

part1 = 2*math.pi*(R*R)

part2 = (math.sin(lat1)) - (math.sin(lat2))

part3 = (lon1 - lon2) / 360

answer = part1 * part2 * part3
``````

Giving an answer of ~456 million ^2. The problem is, the earth's area is ~510km^2. So the result is over 10% out.

• Do you really need a projection here? If all you have is bounding boxes (spherical rectangles), there should be a fairly simple formula based on the lower left and upper right corners. – barrycarter Dec 19 '16 at 17:10
• @barrycarter - that would probably fit the bill, but I have no idea where to begin. The maths stack exchange has a wall of gibberish as an answer to that question - math.stackexchange.com/questions/1205927/… - unless there's a nice pre-written python implementation of that, the simplest solution from an implementation standpoint seems to be to use standard libraries to reproject and calculate, even if slightly overkill (I don't need the ellipsoid level accuracy after all). – GIS-Jonathan Dec 19 '16 at 17:15
• mathforum.org/library/drmath/view/63767.html has a cleaner explanation. – barrycarter Dec 19 '16 at 17:20
• Absolute value. – barrycarter Dec 19 '16 at 17:47
• The problem with bounding boxes is that they're likely sparse--polygon only contains 4 (or 5) points. For larger ones, projecting them into an equal area projection means the output box won't be accurate. If you can densify them, then project them...Densifying can also be problematic, densify along latitude lines, or along geodesic lines? – mkennedy Dec 19 '16 at 20:26

• Thanks. I've tested this out with Gall-Peters because it was the only cylindrical equal-area I could find on spatialreference.org. It produces the correct result for the entire globe, and one degree^2 at the equator. For anyone curious, the pyproj.Proj is: `pyproj.Proj("+proj=cea +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")` – GIS-Jonathan Dec 20 '16 at 17:09