Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
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. –  nmpeterson Apr 14 at 16:23

3 Answers 3

up vote 6 down vote accepted

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)
share|improve this answer
    
Alright thanks a lot for your help! –  Below the Radar Apr 14 at 16:27
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. –  blah238 Apr 14 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. –  blord-castillo Apr 14 at 16:44
    
Isn't Destroy called automatically when Python garbage collects the Python object? –  jpmc26 Apr 14 at 19:52
    
I think it should be, but I noticed in ogr docs that Destroy is normally explicitly called. –  blord-castillo Apr 14 at 20:19

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)
share|improve this answer
    
Is shapely faster than the ogr solution of blord-castillo? –  Below the Radar Apr 14 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() –  blord-castillo Apr 14 at 16:50

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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