def plot_station_mappings(mapping_results):
    """ Plot a list of mapping results on a map.

    Requires matplotlib and cartopy.

    Parameters
    ----------
    mapping_results : list of MappingResult objects
        Mapping results to plot
    """
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        raise ImportError('Plotting requires matplotlib.')
    try:
        import cartopy.crs as ccrs
        import cartopy.feature as cfeature
    except ImportError:
        raise ImportError('Plotting requires cartopy.')
    lats = []
    lngs = []
    t_lats = []
    t_lngs = []
    n_discards = 0
    for mapping_result in mapping_results:
        if not mapping_result.is_empty():
            lat, lng = mapping_result.isd_station.coords
            t_lat, t_lng = map(float, mapping_result.target_coords)
            lats.append(lat)
            lngs.append(lng)
            t_lats.append(t_lat)
            t_lngs.append(t_lng)
        else:
            n_discards += 1
    print('Discarded {} empty mappings'.format(n_discards))
    fig = plt.figure(figsize=(60, 60))
    ax = plt.subplot(1, 1, 1, projection=ccrs.Mercator())
    all_lngs = lngs + t_lngs
    all_lats = lats + t_lats
    x_max = max(all_lngs)
    x_min = min(all_lngs)
    x_diff = x_max - x_min
    y_max = max(all_lats)
    y_min = min(all_lats)
    y_diff = y_max - y_min
    x_pad = 0.1 * x_diff
    y_pad = 0.1 * y_diff
    left = x_min - x_pad
    right = x_max + x_pad
    bottom = y_min - y_pad
    top = y_max + y_pad
    width_ratio = 2.0
    height_ratio = 1.0
    if (right - left) / (top - bottom) > height_ratio / width_ratio:
        goal = (right - left) * height_ratio / width_ratio
        diff = goal - (top - bottom)
        bottom = bottom - diff / 2.0
        top = top + diff / 2.0
    else:
        goal = (top - bottom) * width_ratio / height_ratio
        diff = goal - (right - left)
        left = left - diff / 2.0
        right = right + diff / 2.0
    left = max(left, -179.9)
    right = min(right, 179.9)
    bottom = max([bottom, -89.9])
    top = min([top, 89.9])
    ax.set_extent([left, right, bottom, top])
    ax.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '50m',
        edgecolor='face', facecolor=cfeature.COLORS['water']))
    ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m',
        edgecolor='face', facecolor=cfeature.COLORS['land']))
    ax.add_feature(cfeature.NaturalEarthFeature('cultural',
        'admin_0_boundary_lines_land', '50m', edgecolor='black', facecolor=
        'none'))
    ax.add_feature(cfeature.NaturalEarthFeature('physical', 'lakes', '50m',
        edgecolor='face', facecolor=cfeature.COLORS['water']))
    ax.add_feature(cfeature.NaturalEarthFeature('physical', 'coastline',
        '50m', edgecolor='black', facecolor='none'))
    ax.plot([lngs, t_lngs], [lats, t_lats], color='k', linestyle='-',
        transform=ccrs.Geodetic(), linewidth=0.3)
    ax.plot(lngs, lats, 'bo', markersize=1, transform=ccrs.Geodetic())
    plt.title('Location to weather station mapping')
    plt.show()
