python4oceanographers

Turning ripples into waves

Better colormaps with brewer2mpl

We already used the brewer2mpl module in previous posts. However, we did not showed it in detail. This post is just to draw attention to this interesting module that can access the Colorbrewer2.0 website and convert colormaps directly to matplotlib format.

Let's start importing by cartopy, iris, pyplot, and numpy.

In [2]:
import numpy as np
import matplotlib.pyplot as plt

import iris
import iris.plot as iplt
import cartopy.crs as ccrs
import iris.quickplot as qplt

from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

We will plot a topography map using the etopo1 data from the NOAA opendap server. With iris we can slice the data before downloading the it. Finally we will "print" our cube to check the metadata.

In [3]:
def get_cube(url):
    coord_values = {'latitude':lambda cell: -25 <= cell <= -22,
                    'longitude': lambda cell: -46 <= cell <= -42}
    constraint = iris.Constraint(coord_values=coord_values)
    bathy = iris.load_cube(url, constraint)
    return bathy

try:
    bathy = get_cube('/home/filipe/00-NOBKP/OcFisData/ETOPO1_Bed_g_gmt4.grd')
except:
    bathy = get_cube('http://www.ngdc.noaa.gov/thredds/dodsC/relief/ETOPO1/thredds/ETOPO1_Bed_g_gmt4.nc')


print(bathy)
writing [iris.fileformats._pyke_rules.compiled_krb]/fc_rules_cf_fc.py
writing [iris.fileformats._pyke_rules.compiled_krb]/compiled_pyke_files.py
/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/iris/fileformats/cf.py:1139: UserWarning: NetCDF default loading behaviour currently does not expose variables which define reference surfaces for dimensionless vertical coordinates as independent Cubes. This behaviour is deprecated in favour of automatic promotion to Cubes. To switch to the new behaviour, set iris.FUTURE.netcdf_promote to True.
  warnings.warn(msg)

z / (meters)                        (latitude: 181; longitude: 241)
     Dimension coordinates:
          latitude                           x               -
          longitude                          -               x
     Attributes:
          Conventions: CF-1.4
          GMT_version: 4.5.1 [64-bit]
          actual_range: [-10898.   8271.]
          geospatial_lat_max: 90.0
          geospatial_lat_min: -90.0
          geospatial_lat_resolution: 0.016667
          geospatial_lat_units: degrees_north
          geospatial_lon_max: 180.0
          geospatial_lon_min: -180.0
          geospatial_lon_resolution: 0.016667
          geospatial_lon_units: degrees_east
          history: grdreformat ETOPO1_Bed_g_gmt4.grd ETOPO1_Bed_g_gmt4_coards.grd
          node_offset: 0
          positive: up
          title: ETOPO1_Bed_g.int.grd

Now we will use brewer2mpl to create two color pallets, one with Greens for land, and another with Blues for the ocean. We reversed the ocean pallet because we want the darker blue at deeper depths. In addition, we will cap the first two green levels to get a darker shade of green at the highest heights. Finally, we concatenate both colors and convert them to matplotlib values.

In [4]:
from brewer2mpl import brewer2mpl

land = brewer2mpl.get_map('Greens', 'sequential', 9)
ocean = brewer2mpl.get_map('Blues', 'sequential', 7, reverse=True)

colors = np.array(ocean.mpl_colors + land.mpl_colors[2:])

levels = [-4000, -2500, -1000, -700, -400, -145, -10, 0, 10, 145, 400, 800, 1200, 1600]


fig, ax = plt.subplots(figsize=(6, 6),
                       subplot_kw=dict(projection=ccrs.PlateCarree()))
qplt.contourf(bathy, levels, colors=colors, extend='both')
ax.coastlines('10m', color='k')
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                    linewidth=1.5, color='gray', alpha=0.5, linestyle='--')
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
iplt.citation('Inspired by\nhttp://nbviewer.ipython.org/5254237')
/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/matplotlib/artist.py:221: MatplotlibDeprecationWarning: This has been deprecated in mpl 1.5, please use the
axes property.  A removal date has not been set.
  warnings.warn(_get_axes_msg, mplDeprecation, stacklevel=1)

Much easier (and more readable) then the old MatlabTM way which still needs to perform an overlay of the subplots making the code look even worse than with just the clim trick!

In [5]:
HTML(html)
Out[5]:

This post was written as an IPython notebook. It is available for download or as a static html.

Creative Commons License
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

Comments