After generating some random data, it was obvious that the bounds that I chose did not work with this projection (red lines). Using map.drawgreatcircle(), I first visualized where I wanted the bounds while zoomed over the projection of random data.
I corrected the longitude by using the longitudinal difference at the southern most latitude (blue horizontal line).
I determined the latitudinal range using the Pythagorean theorem to solve for the vertical distance, knowing the distance between the northern most longitudinal bounds, and the central southernmost point (blue triangle).
def centerMap(lats,lons,scale):
#Assumes -90 < Lat < 90 and -180 < Lon < 180, and
# latitude and logitude are in decimal degrees
earthRadius = 6378100.0 #earth's radius in meters
northLat = max(lats)
southLat = min(lats)
westLon = max(lons)
eastLon = min(lons)
# average between max and min longitude
lon0 = ((westLon-eastLon)/2.0)+eastLon
# a = the height of the map
b = sd.spheredist(northLat,westLon,northLat,eastLon)*earthRadius/2
c = sd.spheredist(northLat,westLon,southLat,lon0)*earthRadius
# use pythagorean theorom to determine height of plot
mapH = pow(pow(c,2)-pow(b,2),1./2)
arcCenter = (mapH/2)/earthRadius
lat0 = sd.secondlat(southLat,arcCenter)
# distance between max E and W longitude at most souther latitude
mapW = sd.spheredist(southLat,westLon,southLat,eastLon)*earthRadius
return lat0,lon0,mapW*scale,mapH*scale
lat0center,lon0center,mapWidth,mapHeight = centerMap(dataLat,dataLon,1.1)
The lat0 (or latitudinal center) in this case is therefore the point half-way up the height of this triangle, which I solved using John Cooks method, but for solving for an unknown coordinate while knowing the first coordinate (the median longitude at the southern boundary) and the arc length (half that of the total height).
def secondlat(lat1, arc):
degrees_to_radians = math.pi/180.0
lat2 = (arc-((90-lat1)*degrees_to_radians))*(1./degrees_to_radians)+90
return lat2
Update: The above function, as well as the distance between two coordinates can be achieved with higher accuracy using the pyproj Geod class methods geod.fwd() and geod.inv(). I found this in Erik Westra's Python for Geospatial Development, which is an excellent resource.
Update: I have now verified that this also works for Lambert Conformal Conic (lcc) projections.