docs: 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