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?
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 print(answer)
Giving an answer of ~456 million ^2. The problem is, the earth's area is ~510km^2. So the result is over 10% out.