2

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
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.

Improve this question
8
  • 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.
    – user1462
    Dec 19, 2016 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, 2016 at 17:15
  • mathforum.org/library/drmath/view/63767.html has a cleaner explanation.
    – user1462
    Dec 19, 2016 at 17:20
  • Absolute value.
    – user1462
    Dec 19, 2016 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, 2016 at 20:26

1 Answer 1

Reset to default
1

Cylindrical eqal area should work quite well because the parallels remain parallel (good for bottom and top of bounding boxes) and the dateline is a single vertical line (not so difficult to split).

Improve this answer
1
  • 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, 2016 at 17:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.