Minimum bounding box for set of bounding boxes (longitude calculation)

boundingbox.py

import numpy as np
import numpy.ma as ma

# this is an implementation of:
# http://gis.stackexchange.com/questions/17788/how-to-compute-the-bounding-box-of-multiple-layers-in-lat-long/17987#17987

def minimumBoundingBox(lons):
    
    xs = np.sort(lons.ravel())
    xs = np.concatenate((xs,[xs[0]+360]))
    
    lonsUnwrapped = np.rad2deg(np.unwrap(np.deg2rad(lons)))
    
    coversInterval = np.zeros(len(xs)-1, dtype=bool)
    for i in xrange(1, len(xs)):
        for bb in lonsUnwrapped:
            if bb[0] <= xs[i-1] and bb[1] >= xs[i]:
                coversInterval[i-1] = True
                break
    
    intervalLengths = xs[1:] - xs[:-1]
    gapLengths = ma.masked_array(intervalLengths, coversInterval)
    biggestGapIdx = np.argmax(gapLengths) # ignoring other possible maximums
    
    lonWest, lonEast = xs[biggestGapIdx+1], xs[biggestGapIdx]
    # wrap into [-180,180]
    lonWest = np.rad2deg((np.deg2rad(lonWest) + np.pi) % (2 * np.pi ) - np.pi)
    lonEast = np.rad2deg((np.deg2rad(lonEast) + np.pi) % (2 * np.pi ) - np.pi)
    
    print lons
    print xs
    print lonsUnwrapped 
    print coversInterval
    print gapLengths
    print biggestGapIdx
    
    return lonWest, lonEast
    
if __name__ == '__main__':
    #lons = np.array([[-81,-77],[77,156],[-149,-90],[-69,-36]])
    lons = np.array([[170,-170],[160,175]])
    print minimumBoundingBox(lons)