def plot_mag_map(fignum, element, lons, lats, element_type, cmap='coolwarm', lon_0=0, date="", contours=False, proj='PlateCarree'):
    """
    makes a color contour map of geomagnetic field element

    Parameters
    ____________
    fignum : matplotlib figure number
    element : field element array from pmag.do_mag_map for plotting
    lons : longitude array from pmag.do_mag_map for plotting
    lats : latitude array from pmag.do_mag_map for plotting
    element_type : [B,Br,I,D] geomagnetic element type
        B : field intensity
        Br : radial field intensity
        I : inclinations
        D : declinations
    Optional
    _________
    contours : plot the contour lines on top of the heat map if True
    proj : cartopy projection ['PlateCarree','Mollweide']
           NB: The Mollweide projection can only be reliably with cartopy=0.17.0; otherwise use lon_0=0.  Also, for declinations, PlateCarree is recommended.
    cmap : matplotlib color map - see https://matplotlib.org/examples/color/colormaps_reference.html for options
    lon_0 : central longitude of the Mollweide projection
    date : date used for field evaluation,
           if custom ghfile was used, supply filename

    Effects
    ______________
    plots a color contour map with  the desired field element
    """

    if not has_cartopy:
        print('This function requires installation of cartopy')
        return
    from matplotlib import cm
    if lon_0 == 180:
        lon_0 = 179.99
    if lon_0 > 180:
        lon_0 = lon_0-360.
    lincr = 1
    if type(date) != str:
        date = str(date)
    if proj == 'PlateCarree':
        fig = plt.figure(fignum)
        ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=lon_0))
    if proj == 'Mollweide':
        fig = plt.figure(fignum)
        # this issue is fixed in >=0.17
        if not LooseVersion(Cartopy.__version__) > LooseVersion('0.16.0'):
            if lon_0 != 0:
                print('This projection requires lon_0=0')
                return

        ax = plt.axes(projection=ccrs.Mollweide(central_longitude=lon_0))
    xx, yy = np.meshgrid(lons, lats)
    levmax = 5*round(element.max()/5)+5
    levmin = 5*round(element.min()/5)-5
    if element_type == 'Br' or element_type == 'B':
        plt.contourf(xx, yy, element,
                     levels=np.arange(levmin, levmax, 1),
                     cmap=cmap, transform=ccrs.PlateCarree())
        cbar = plt.colorbar(orientation='horizontal')
        if contours:
            plt.contour(xx, yy, element, levels=np.arange(levmin, levmax, 10),
                        colors='black', transform=ccrs.PlateCarree())

        if element_type == 'Br':
            plt.title('Radial field strength ($\mu$T): '+date)
        else:
            plt.title('Total field strength ($\mu$T): '+date)
    if element_type == 'I':
        plt.contourf(xx, yy, element,
                     levels=np.arange(-90, 90, lincr),
                     cmap=cmap, transform=ccrs.PlateCarree())
        cbar = plt.colorbar(orientation='horizontal')
        if contours:
            plt.contour(xx, yy, element, levels=np.arange(-80, 90, 10),
                        colors='black', transform=ccrs.PlateCarree())
        plt.title('Field inclination: '+date)
    if element_type == 'D':
        plt.contourf(xx, yy, element,
                     levels=np.arange(-180, 180, 10),
                     cmap=cmap, transform=ccrs.PlateCarree())
        cbar = plt.colorbar(orientation='horizontal')
        if contours:
            plt.contour(xx, yy, element, levels=np.arange(-180, 180, 10),
                        colors='black', transform=ccrs.PlateCarree())
        plt.title('Field declination: '+date)
    ax.coastlines()
    ax.set_global()
    return ax