[GIS] Forcing Lat/Lon Bounding Box to appear on LCC Projection Map in matplotlib-Basemap

cartographymatplotlib

I'm trying to do some plots of data on a basemap of Australia using
matplotlib-Basemap. I'm finding the interface for creating the basemap
to be tricky. Specifically, I have lat/lon definitions for the Australian
region (including Tasmania), but the map set up using the Basemap API excludes
parts of the desired region. Here's a small sample code, which invokes
matplotlib-Basemap's Basemap with the BoM's LCC projection parameters. If
you run this, you'll find the northern portion of the Cape York Peninsula and
the western portion of WA are missing. I'm hoping to do this automatically for
Australia and for two of its states for which the standard map projection is
also LCC (Tasmania and Victoria).

How can I get a basemap that will have as a guarantee that it contains all of a
lat/lon region defined by a lat/lon bounding box?

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

"""
Equatorial/Polar minor radii for WGS84 ellipsoid
approximation of the geoid
"""
rEquat = 6378137.00
rPolar = 6356752.3142
"""
Standard latitudes for LCC projection.  Used in this
projection by Australia's BoM.
"""
trueLat1 = -10.
trueLat2 = -40.

"""
lower-left/upper-right corners for the Australian domain.
"""
ozMinLat = -43.575
ozMinLon = 112.925
ozMaxLat = -10.075
ozMaxLon = 153.575

"""
Central lat/lon coordinates.
"""
centerLat = 0.5 * (ozMinLat + ozMaxLat)
centerLon = 0.5 * (ozMinLon + ozMaxLon)

"""
CASE 1:
BoM's projection parameters.  Cuts off W. edge of 
WA and N end of Cape York Peninsula.
"""
m = Basemap(llcrnrlon=ozMinLon, llcrnrlat=ozMinLat, 
            urcrnrlon=ozMaxLon, urcrnrlat= ozMaxLat,
            rsphere=(rEquat,rPolar), resolution='c',
            area_thresh=1000.,projection='lcc',
            lat_1=trueLat1, lat_2=trueLat2, 
            lat_0=centerLat, lon_0=centerLon)

m.drawcoastlines()            
plt.show()

Best Answer

As you note in the comments in your code, the max/min latitudes and longitudes don't define a bounding box, but the geodetic lat/long of the lower left and upper right corner of the image. As you can see from the plot below, Basemap completely respects these boundaries (plotted for longitude only). However, as parallels and meridians aren't parallel to each other when projected in the Lambert Conformal projection, it is impossible to observe the same values for the top left corner as well, or indeed along the left or top edge.

This is with the following Basemap call:

m = Basemap(llcrnrlon=ozMinLon, llcrnrlat=ozMinLat, 
        urcrnrlon=ozMaxLon, urcrnrlat= ozMaxLat,
        rsphere=(rEquat,rPolar), resolution='c',
        area_thresh=1000.,projection='lcc', lat_0=centerLat, lon_0=centerLon)

Test run of user code with meridians plotted

Note that the Basemap documentation indicates:

For most map projections, the map projection region can either be specified by setting these keywords:

llcrnrlon longitude of lower left hand corner of the desired map domain (degrees).
llcrnrlat latitude of lower left hand corner of the desired map domain (degrees).
urcrnrlon longitude of upper right hand corner of the desired map domain (degrees).
urcrnrlat latitude of upper right hand corner of the desired map domain (degrees).
or these

width width of desired map domain in projection coordinates (meters).
height height of desired map domain in projection coordinates (meters).
lon_0 center of desired map domain (in degrees).
lat_0 center of desired map domain (in degrees).

So one possibility would be to specify the center and pass in the total width and height (NS and EW extent) of the Australian landmass.