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

<...>

def main():
# The rest of the main program
<...>

    # Plot data on a map
    if maps:
        print('\nPlotting {:d} sample location maps'.format(tot_maps))
        # plot maps for MI and derived parameters
        if (majorions) :
            tot_maps += MI_maps
            for analyte in range(MI_maps) :
                map_data = {'param' : measurements_MI[MI_map_list[analyte]],
                            'lats' : samp_info_MI['Lat_St'],
                            'lons' : samp_info_MI['Lon_St'],
                            'projection' : 'globe-cyl',
                            'plot_lim' : [0.,0.,0.,0.],
                            'fig_size' : gl_cyl_map_size,
                            'file_form' : picformat}

                plot_data(map_data, MI_map_list[analyte], MI_title_list[analyte])


#===============================================================================
def plot_data(map_data, figname, title):
    """
    Create sample position maps for each of the measured parameters

    Arguments
    ---------
    map_data : dict
        dictionary with plot parameters
    figname : string
        Name of the saved figure
    title : string
        Figure title
    """
    datain = map_data['param']
    plats = map_data['lats']
    plons = map_data['lons']
    proj = map_data['projection']
    box = map_data['plot_lim']
    psize = map_data['fig_size']
    picformat = map_data['file_form']

    idx = ~np.isnan(datain)
    longitudes = plons[idx]
    latitudes = plats[idx]

    fig = plt.figure(figsize=psize, dpi=300)
    if proj == 'orthographic':
        gs = GridSpec(100, 100, bottom=0.05, top=0.95, left=0.10, right=0.90)
        ax = fig.add_subplot(gs[20:99,0:99])
    else :
        gs = GridSpec(100, 100, bottom=0.05, top=0.95, left=0.10, right=0.95)
        ax = fig.add_subplot(gs[10:99,0:99])

    if proj == 'globe-cyl':
        s_map = Basemap(projection='cyl', lat_0=0, lon_0=0, resolution = 'l')
    elif proj == 'cylinder':
        s_map = Basemap(projection='cyl', llcrnrlat=box[0], urcrnrlat=box[1],
                        llcrnrlon=box[2], urcrnrlon=box[3], resolution = 'i')
    elif (proj == 'orthographic'):
        s_map = Basemap(projection='ortho', lat_0 = box[0], lon_0 = box[1],
                        resolution = 'l', ax=ax)
    s_map.drawcoastlines()
    s_map.fillcontinents('lightgreen')
    if proj == 'globe-cyl':
        s_map.drawparallels(np.arange(-90., 91., 30.), labels=[1,0,0,0])
        s_map.drawmeridians(np.arange(-180., 181., 60.), labels=[0,0,0,1])
    elif proj == 'cylinder' :
        s_map.drawparallels(np.arange(box[0], box[1]+1, 20.), labels=[1,0,0,0])
        s_map.drawmeridians(np.arange(box[2], box[3]+1, 20.), labels=[0,0,0,1])
    elif proj == 'orthographic':
        s_map.drawparallels(np.arange(-90., 91., 30.))
        s_map.drawmeridians(np.arange(-180., 181., 30.))

    s_map.plot(plons, plats, 'yo', ms=2, latlon='True')
    s_map.plot(longitudes, latitudes, 'ro', ms=6, latlon='True')
    ttl = 'Places where {title} was measured \n(n = {count:g})\n'
    ax.set_title(ttl.format(title=title,
                            count=np.count_nonzero(~np.isnan(longitudes))))

    name = '{}.{}'.format(figname, picformat)
    fig.savefig(name)
    plt.close(fig)

#===============================================================================
# Make sure the main program is executed
#===============================================================================
if __name__ == '__main__':
    main()