3

I am trying to find a way (simplest is the best) to retrieve the largest extent in a list of tuple extents in python.

The list is looking like that:

extents = [(-180.0, -90.0, 180.0, 83.624), (-124.731, 24.956, -66.97, 49.372), (-122.42, - 37.818, 151.207, 52.516)]

Could it be a good idea to sum the absolute value of each coordinates and retrieve the highest?

Or does a geographic library such OGR could be used?

Edit:

I have created Envelope objects with GDAL library:

envelopes = [Envelope(extent) for extent in extents]

Now I have to compare the area of each? Any suggestion?

Can anybody give explanation on this? Thanks a lot.

Improve this question
1
  • 2
    If you're defining extents in degrees, simple math is dangerous: extents of the "same size" (in degrees) are very different near the equator vs near the poles. Commented Apr 14, 2014 at 16:23

3 Answers 3

Reset to default
6

Main problem is getting the area of the extent. I wrote a quick ogr function to do this

def extentArea(extent): #Unpack extent tuple to coordinates minX, minY, maxX, maxY = extent #unpack the tuple #Create empty geometry and add vertices geom = ogr.Geometry(type = ogr.wkbLinearRing) geom.AddPoint_2D(minX,minY) geom.AddPoint_2D(minX,maxY) geom.AddPoint_2D(maxX,maxY) geom.AddPoint_2D(maxX,minY) geom.AddPoint_2D(minX,minY) # Plan to return area, but destroy geometry first area = geom.GetArea() geom.Destroy() return area

Once you have that function, finding the largest extent is simple with max() using extentArea as the sorting key.

>>> extents = [(-180.0, -90.0, 180.0, 83.624), (-124.731, 24.956, -66.97, 49.372), (-122.42, -37.818, 151.207, 52.516)] >>> max(extents, key=extentArea) (-180.0, -90.0, 180.0, 83.624)
Improve this answer
4
  • 2
    I think that nmpeterson's comment still applies here. It would probably be safer to project the extents into an equal-area projection before attempting to compute the area. GetArea() uses Green's theorem which, from my understanding, does not take into account the fact that latitude and longitude are angular, not linear units. Commented Apr 14, 2014 at 16:38
  • Very true. Especially a problem for very large extents like those in the example extents. Although... I am not certain that is an issue when comparing extents if you only need relative size, not exact value. It probably is an issue, just not certain it would be an issue. Commented Apr 14, 2014 at 16:44
  • Isn't Destroy called automatically when Python garbage collects the Python object? Commented Apr 14, 2014 at 19:52
  • I think it should be, but I noticed in ogr docs that Destroy is normally explicitly called. Commented Apr 14, 2014 at 20:19
4

It is easiest with shapely:

from shapely.geometry import box extents = [(-180.0, -90.0, 180.0, 83.624), (-124.731, 24.956, -66.97, 49.372), (-122.42, -37.818, 151.207, 52.516)] for i in extents: a = box(i[0],i[1],i[2],i[3]) print i, a.area (-180.0, -90.0, 180.0, 83.623999999999995) 62504.64 (-124.73099999999999, 24.956, -66.969999999999999, 49.372) 1410.292576 (-122.42, -37.817999999999998, 151.20699999999999, 52.515999999999998) 24717.821418

New

with the solution of blord-castillo:

def extentArea(i): geom = box(i[0],i[1],i[2],i[3]) return geom.area print max(extents, key=extentArea) (-180.0, -90.0, 180.0, 83.623999999999995)
Improve this answer
2
  • Is shapely faster than the ogr solution of blord-castillo? Commented Apr 14, 2014 at 16:47
  • Both should be pretty fast. Shapely might be faster because of the box case. You would need a large number of operations to notice a difference. (And if you are doing large numbers of operations, the reprojection problem is bigger) It does introduce an additional dependency though. Once again, you just substitute for the getArea function with a shapely based one, and then use that function as your key on max() Commented Apr 14, 2014 at 16:50
3

If you want to do this "right" (taking into account the fact that latitude and longitude are angular units, and using an ellipsoid as a model of the earth's shape), you can try using the geographiclib library, which is a Python version of Charles Karney's Algorithms for geodesics.

See also the Wikipedia page Geodesics on an ellipsoid for a look into some of the math involved.

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.