# make plots of etopo bathymetry/topography data on
# various map projections, drawing coastlines, state and
# country boundaries, filling continents and drawing
# parallels/meridians

# illustrates special-case polar-centric projections.

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

# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
etopo = np.loadtxt('etopo20data.gz')
lons = np.loadtxt('etopo20lons.gz')
lats = np.loadtxt('etopo20lats.gz')

print 'min/max etopo20 data:'
print etopo.min(),etopo.max()

# these are the 4 polar projections
projs = ['laea','stere','aeqd','ortho'] # short names
# long names
projnames = ['Lambert Azimuthal Equal Area','Stereographic','Azimuthal Equidistant','Orthographic']
# loop over hemispheres, make a 4-panel plot for each hemisphere
# showing all four polar projections.
for hem in ['North','South']:
    if hem == 'South':
        lon_0 = 130.
        lon_0_ortho = lon_0 - 180.
        lat_0 = -90.
        bounding_lat = -20.
    elif hem == 'North':
        lon_0 = -90.
        lon_0_ortho = lon_0
        lat_0 = 90.
        bounding_lat = 20.
    # loop over projections, one for each panel of the figure.
    fig = plt.figure(figsize=(8,8))
    npanel = 0
    for proj,projname in zip(projs,projnames):
        npanel = npanel + 1
        if hem == 'South':
            projection = 'sp'+proj
        elif hem == 'North':
            projection = 'np'+proj
        # setup map projection
        # centered on Australia (for SH) or US (for NH).
        if proj == 'ortho':
           m = Basemap(projection='ortho',
                       resolution='c',area_thresh=10000.,lat_0=lat_0,lon_0=lon_0_ortho)
        else:
           m = Basemap(boundinglat=bounding_lat,lon_0=lon_0,\
                       resolution='c',area_thresh=10000.,projection=projection)
        # compute native map projection coordinates for lat/lon grid.
        x,y = m(*np.meshgrid(lons,lats))
        ax = fig.add_subplot(2,2,npanel)
        # make filled contour plot.
        cs = m.contourf(x,y,etopo,20,cmap=plt.cm.jet)
        # draw coastlines.
        m.drawcoastlines()
        # draw parallels and meridians.
        m.drawparallels(np.arange(-80.,90,20.))
        m.drawmeridians(np.arange(0.,360.,60.))
        # draw boundary around map region.
        m.drawmapboundary()
        # draw title.
        plt.title(hem+' Polar '+projname,y=1.05,fontsize=12)
        print 'plotting '+hem+' Polar '+projname+' basemap ...'
plt.show()