2

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()
3

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.

|improve this answer|||||
0

I have a inelegant kludge that works for me ... which I came to by trial-and-error ...

OZMINLAT = -44+1     # (South East Cape Tasmania 43° 38' 40" S 146° 49' 30" E)
OZMAXLAT = -9+1.5    # (Boigu Island 09° 16' S 142° 13' E)
OZMINLON = 112.5-4.1 # (Steep point 26° 09' 5" S 113° 09' 18" E
OZMAXLON = 154-2.5   # (Cape Byron 28° 38' 15" S 153° 38' 14" E)
OZMIDLAT = -25.6
OZMIDLON = 134.35

m = Basemap(llcrnrlon=OZMINLON, llcrnrlat=OZMINLAT, 
    urcrnrlon=OZMAXLON, urcrnrlat=OZMAXLAT, resolution='i',
    lat_0=OZMIDLAT, lon_0=OZMIDLON, 
    area_thresh=100.,projection='lcc')

Resultant map of Australia

|improve this answer|||||

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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