diff --git a/contour_example_gl.py b/contour_example_gl.py
new file mode 100755
index 0000000..6653f7a
--- /dev/null
+++ b/contour_example_gl.py
@@ -0,0 +1,100 @@
+# -*- coding: utf-8 -*-
+"""
+This example demonstrates the use of GLSurfacePlotItem.
+"""
+
+
+## Add path to library (just for examples; you do not need this)
+#import initExample
+
+from pyqtgraph.Qt import QtCore, QtGui
+import pyqtgraph as pg
+import pyqtgraph.opengl as gl
+import numpy as np
+
+## Create a GL View widget to display data
+app = QtGui.QApplication([])
+w = gl.GLViewWidget()
+w.show()
+w.setWindowTitle('pyqtgraph example: GLSurfacePlot')
+w.setCameraPosition(distance=50)
+
+## Add a grid to the view
+g = gl.GLGridItem()
+g.scale(2,2,1)
+g.setDepthValue(10)  # draw grid after surfaces since they may be translucent
+w.addItem(g)
+
+
+## Simple surface plot example
+## x, y values are not specified, so assumed to be 0:50
+z = pg.gaussianFilter(np.random.normal(size=(50,50)), (1,1))
+p1 = gl.GLSurfacePlotItem(z=z, shader='shaded', color=(0.5, 0.5, 1, 1))
+p1.scale(16./49., 16./49., 1.0)
+p1.translate(-18, 2, 0)
+w.addItem(p1)
+
+
+## Saddle example with x and y specified
+x = np.linspace(-8, 8, 50)
+y = np.linspace(-8, 8, 50)
+z = 0.1 * ((x.reshape(50,1) ** 2) - (y.reshape(1,50) ** 2))
+p2 = gl.GLSurfacePlotItem(x=x, y=y, z=z, shader='normalColor')
+p2.translate(-10,-10,0)
+w.addItem(p2)
+
+
+## Manually specified colors
+z = pg.gaussianFilter(np.random.normal(size=(50,50)), (1,1))
+x = np.linspace(-12, 12, 50)
+y = np.linspace(-12, 12, 50)
+colors = np.ones((50,50,4), dtype=float)
+colors[...,0] = np.clip(np.cos(((x.reshape(50,1) ** 2) + (y.reshape(1,50) ** 2)) ** 0.5), 0, 1)
+colors[...,1] = colors[...,0]
+
+p3 = gl.GLSurfacePlotItem(z=z, colors=colors.reshape(50*50,4), shader='shaded', smooth=False)
+p3.scale(16./49., 16./49., 1.0)
+p3.translate(2, -18, 0)
+w.addItem(p3)
+
+
+
+
+## Animated example
+## compute surface vertex data
+cols = 90
+rows = 100
+x = np.linspace(-8, 8, cols+1).reshape(cols+1,1)
+y = np.linspace(-8, 8, rows+1).reshape(1,rows+1)
+d = (x**2 + y**2) * 0.1
+d2 = d ** 0.5 + 0.1
+
+## precompute height values for all frames
+phi = np.arange(0, np.pi*2, np.pi/20.)
+z = np.sin(d[np.newaxis,...] + phi.reshape(phi.shape[0], 1, 1)) / d2[np.newaxis,...]
+
+
+## create a surface plot, tell it to use the 'heightColor' shader
+## since this does not require normal vectors to render (thus we 
+## can set computeNormals=False to save time when the mesh updates)
+p4 = gl.GLSurfacePlotItem(x=x[:,0], y = y[0,:], shader='heightColor', computeNormals=False, smooth=False)
+p4.shader()['colorMap'] = np.array([0.2, 2, 0.5, 0.2, 1, 1, 0.2, 0, 2])
+p4.translate(10, 10, 0)
+w.addItem(p4)
+
+index = 0
+def update():
+    global p4, z, index
+    index -= 1
+    p4.setData(z=z[index%z.shape[0]])
+    
+timer = QtCore.QTimer()
+timer.timeout.connect(update)
+timer.start(30)
+
+## Start Qt event loop unless running in interactive mode.
+if __name__ == '__main__':
+    import sys
+    if (sys.flags.interactive != 1) or not hasattr(QtCore, 'PYQT_VERSION'):
+        QtGui.QApplication.instance().exec_()
+
diff --git a/contour_sample.py b/contour_sample.py
new file mode 100644
index 0000000..02ddca1
--- /dev/null
+++ b/contour_sample.py
@@ -0,0 +1,42 @@
+from pyqtgraph.Qt import QtGui, QtCore
+import numpy as np
+import pyqtgraph as pg
+import sys
+
+# Setup
+app = QtGui.QApplication([])
+pg.setConfigOption('background', 'w')
+pg.setConfigOption('foreground', 'k')
+
+win = pg.PlotWidget()
+layout = pg.GraphicsLayout()
+win.setCentralItem(layout)
+ax = pg.PlotItem()
+layout.addItem(ax)
+
+# Generate data
+x = np.linspace(10, 16.28, 30)
+y = x[:]
+xx, yy = np.meshgrid(x, y)
+z = np.sin(xx) + np.cos(yy)
+
+# Add data
+ax.setXRange(x.min(), x.max())
+ax.setYRange(y.min(), y.max())
+
+c = pg.IsocurveItem(data=z, level=0.5, pen='r', axisOrder='row-major')
+img = pg.ImageItem(z, axisOrder='row-major')
+img.translate(x.min(), y.min())
+img.scale((x.max() - x.min()) / img.width(), (y.max() - y.min()) / img.height())
+ax.addItem(img)
+
+# c.setParentItem(img)
+# https://stackoverflow.com/a/51109935/6605826
+c.translate(x.min(), y.min())
+c.scale((x.max() - x.min()) / np.shape(z)[0], (y.max() - y.min()) / np.shape(z)[1])
+
+ax.addItem(c)
+
+# Finish up
+win.show()
+sys.exit(app.exec_())
diff --git a/contour_sample2.py b/contour_sample2.py
new file mode 100644
index 0000000..06345fa
--- /dev/null
+++ b/contour_sample2.py
@@ -0,0 +1,32 @@
+from pyqtgraph.Qt import QtGui, QtCore
+import numpy as np
+import pyqtgraph as pg
+import sys
+
+# Setup
+app = QtGui.QApplication([])
+pg.setConfigOption('background', 'w')
+pg.setConfigOption('foreground', 'k')
+
+win = pg.PlotWidget()
+layout = pg.GraphicsLayout()
+win.setCentralItem(layout)
+ax = pg.PlotItem()
+layout.addItem(ax)
+
+# Generate data
+x = np.linspace(0, 6.28, 30)
+y = x[:]
+xx, yy = np.meshgrid(x, y)
+z = np.sin(xx) + np.cos(yy)
+
+# Add data
+#ax.setXRange(x.min(), x.max())
+#ax.setYRange(y.min(), y.max())
+c = pg.IsocurveItem(data=z, level=0.5, pen='r', axisOrder='row-major')
+# c.setParentItem(ax)  # This doesn't work, of course
+ax.addItem(c)
+
+# Finish up
+win.show()
+sys.exit(app.exec_())
diff --git a/pgmpl/contour.py b/pgmpl/contour.py
index fb5ccfa..8124c17 100644
--- a/pgmpl/contour.py
+++ b/pgmpl/contour.py
@@ -45,6 +45,14 @@ def __init__(self, ax, *args, **kwargs):
         self.nchunk = kwargs.pop('nchunk', 0)
         self.x, self.y, self.z, self.levels = self.choose_xyz_levels(*args)
         self.auto_range(self.z)
+
+        # Get the boundary path
+        self.boundary = [
+            (self.x.min(), self.y.min()),
+            (self.x.max(), self.y.min()),
+            (self.x.max(), self.y.max()),
+            (self.x.min(), self.y.max()),
+        ]
         self.draw()
 
         return
@@ -52,7 +60,7 @@ def __init__(self, ax, *args, **kwargs):
     def auto_range(self, z):
         pad = (z.max() - z.min()) * 0.025
         self.vmin = z.min()+pad if self.vmin is None else self.vmin
-        self.vmax = z.max()-pad if self.vmax is None else self.vmax
+        self.vmax = z.max() + (pad if self.filled else -pad) if self.vmax is None else self.vmax
 
     def auto_pick_levels(self, z, nlvl=None):
         """
@@ -62,7 +70,7 @@ def auto_pick_levels(self, z, nlvl=None):
             Number of levels; set to some arbitrary default if None
         :return: array
         """
-        nlvl = 5 if nlvl is None else nlvl
+        nlvl = 8 if nlvl is None else nlvl
         self.auto_range(z)
         return np.linspace(self.vmin, self.vmax, nlvl)
 
@@ -89,11 +97,11 @@ def choose_xyz_levels(self, *args):
         levels = lvlinfo if ((lvlinfo is not None) and np.iterable(lvlinfo)) else self.auto_pick_levels(z, lvlinfo)
 
         if x is None:
-            x, y = np.arange(np.shape(z)[0]), np.arange(np.shape(z)[1])
+            x, y = np.arange(np.shape(z)[1]), np.arange(np.shape(z)[0])
 
-        return x, y, z, levels
+        return x, y, z.T, levels
 
-    def extl(self, v):
+    def _extl(self, v):
         """
         Casts input argument as a list and ensures it is at least as long as levels
         :param v: Some variable
@@ -107,15 +115,235 @@ def draw(self):
             self.colors = color_map_translator(
                 self.levels, **{a: self.__getattribute__(a) for a in ['alpha', 'cmap', 'norm', 'vmin', 'vmax']})
         else:
-            self.colors = self.extl(self.colors)
+            self.colors = self._extl(self.colors)
 
         if self.filled:
-            raise NotImplementedError('Filled contours with contourf not yet finished.')
+            self.draw_filled()
         else:
             self.draw_unfilled()
 
+    def _scale_contour_lines(self, lines):
+        """
+        Translates and stretches contour lines
+        :param lines: A list of path segments. e.g. [[(x, y), (x, y)], [(x, y), (x, y)]]
+        :return: A list of path segments, shifted and stretched to fit the x, y data range.
+        """
+        x0, y0, x1, y1 = self.x.min(), self.y.min(), self.x.max(), self.y.max()
+        for i in range(len(lines)):
+            for j in range(len(lines[i])):
+                newx = lines[i][j][0] * (x1 - x0) / np.shape(self.z)[0] + x0
+                newy = lines[i][j][1] * (y1 - y0) / np.shape(self.z)[1] + y0
+                lines[i][j] = (newx, newy)
+        return lines
+
+    def _detect_direction(self, x, y):
+        """
+        Tries to decide whether path goes CCW or not
+
+        :param x: array of numbers
+        :param y: array of numbers
+        :return: True if path appears to run CCW
+        """
+
+        if len(x) < 2:
+            return True  # No need to reverse if 0 or 1 element
+
+        x0, y0 = np.mean([self.x.min(), self.x.max()]), np.mean([self.y.min(), self.y.max()])
+        theta = np.arctan2(y-y0, x-x0)
+
+        # Find theta increment
+        dth = np.diff(theta)
+
+        # Unwrap theta increment so there are no jumps as theta loops around
+        mdth = np.median(dth)
+        jumps = (np.sign(dth) == -np.sign(mdth)) & (abs(dth) > (2*abs(mdth)))
+        for jump in np.where(jumps)[0]:
+            dth[jump] -= 2*np.pi * np.sign(dth[jump])
+
+        return np.mean(dth) > 0
+
+    def _find_which_edge(self, xy_point):
+        print(xy_point, 'xy_point')
+        a = np.where([xy_point[0] == self.x.max(), xy_point[1] == self.y.max(),
+                      xy_point[0] == self.x.min(), xy_point[1] == self.y.min()])[0]
+        if len(a):
+            return a[0]
+        else:
+            return 4
+
+    def _join_lines(self, lines):
+        """
+        Joins segments of a broken path. Use for contours which cross out of bounds and back in. Has to find the right
+        direction of path segments and join them in the right order.
+
+        :param lines: List of path segments. e.g. [[(x, y), (x,y)], [(x, y), (x, y)]]
+        :return: Flattened list of correctly joined path segments e.g. [(x, y), (x,y), (x, y), (x, y)], running CCW
+        """
+        for i in range(len(lines)):
+            # Make sure all segments run CCW within themselves
+            x, y = map(list, zip(*lines[i]))
+            if not self._detect_direction(x, y):
+                lines[i] = lines[i][::-1]
+        # Make sure the set of start points of each segment runs CCW
+        x1 = np.array([line[0][0] for line in lines])
+        y1 = np.array([line[0][1] for line in lines])
+        if not self._detect_direction(x1, y1):
+            lines = lines[::-1]
+
+        if not lines:
+            print('return early because not lines:', lines)
+            return tolist(np.array(self.boundary)[np.array([0, 1, 2, 3, 0])])
+
+        corners = np.roll(self.boundary, -2)
+
+        def get_more_corners(next_e0_, edge_):
+            print('edge: {}, next_edge start: {}'.format(edge_, next_e0_))
+            corners_ = []
+            if (next_e0_ - edge_) >= 1:
+                print('new corner = ', corners[edge_])
+                corners_ += [tuple(corners[edge_]), tuple(corners[edge_+1])]
+                print('Add corner {}'.format(corners[edge_], corners[edge_+1]))
+            if (next_e0_ - edge_) >= 2:
+                corners_ += [tuple(corners[edge_+2])]
+                print('Add corner {}'.format(corners[edge_+2]))
+            if (next_e0_ - edge_) >= 3:
+                print('Add corner {}'.format(corners[edge_+3]))
+                corners_ += [tuple(corners[edge_+3])]#, tuple(corners[edge_])]
+            return corners_
+
+        oneline = lines.pop(0)
+        #if not lines:
+        #    return oneline
+
+        while len(lines):
+            # Find which edge each segment starts/ends on
+            edge = self._find_which_edge(oneline[-1])
+            edge0 = np.array([self._find_which_edge(line[0]) for line in lines])
+            if not any(edge0 >= edge):
+                edge -= 4
+            eligible = [line for line, e0 in zip(lines, edge0) if e0 >= edge]
+            next_line = eligible[edge0[edge0 >= edge].argmin()]
+            next_e0 = edge0[edge0 >= edge].min()
+            oneline += get_more_corners(next_e0, edge)
+
+            lines.remove(next_line)
+            oneline += next_line
+
+        e0 = self._find_which_edge(oneline[0])
+        e1 = self._find_which_edge(oneline[-1])
+        if e1 > e0:
+            e1 -= 4
+        print('e0, e1', e0, e1)
+        more_corners = get_more_corners(e0, e1)
+        print('more corners = ', more_corners)
+        oneline += get_more_corners(e0, e1)
+        oneline += [oneline[0]]
+        return oneline
+
+    def _close_curve(self, line):  # Remove during cleanup -----------------------------------
+        """
+        OLD IDEA; DON'T USE ANYMORE
+        Closes a curve which may be open between the two end points because it intersects the edge of the plot
+        :param line: List of vertices along a CCW path. e.g. [(x, y), (x, y), ...] such as output by _join_lines()
+        :return: List of vertices along a closed CCW path
+        """
+
+        if not line:
+            # Empty line; nothing to do
+            return self.boundary
+
+        if all(np.atleast_1d(line[0] == line[-1])):
+            # Path is already closed; return it with no changes
+            return line
+
+        # Determine which endpoints are on which edges
+        edge0 = np.append(line[0][0] == np.array([self.x.min(), self.x.max()]),
+                          line[0][1] == np.array([self.y.min(), self.y.max()]))
+        edge1 = np.append(line[-1][0] == np.array([self.x.min(), self.x.max()]),
+                          line[-1][1] == np.array([self.y.min(), self.y.max()]))
+        same_edge = edge0 == edge1
+
+        if (not any(edge0) and not any(edge1)) or all(same_edge):
+            # Endpoints are not on the edges, or are on the same edge. Just close the loop.
+            newline = line + [line[0]]
+        else:
+            # Endpoints are not on the same edge; complicated closure
+
+            # Find which boundary path points are between the endpoints of the curve
+            x0, y0 = np.mean([self.x.min(), self.x.max()]), np.mean([self.y.min(), self.y.max()])
+            theta0 = np.arctan2(line[0][1] - y0, line[0][0] - x0)
+            theta1 = np.arctan2(line[-1][1] - y0, line[-1][0] - x0)
+            thetab = np.array([np.arctan2(b[1] - y0, b[0] - x0) for b in self.boundary])
+            # Make sure the thing wraps the right way; we are continuing from point -1 back to point 0
+            theta0 = theta0 + 2*np.pi if theta0 < theta1 else theta0
+            thetab = np.array([b + (2*np.pi if b < theta1 else 0) for b in thetab])
+            # Add in corners of the data range boundary to complete the curve
+            newline = line + tolist(np.array(self.boundary)[thetab < theta0])
+            # And finally, close it
+            newline += [line[0]]
+        return newline
+
+    def draw_filled(self):
+        pens = [setup_pen_kw(penkw=dict(color=self.colors[i])) for i in range(len(self.levels))]
+        use_pen = pg.mkPen(color='k')
+        curves = [None] * (len(self.levels))
+        joined_lines = [None] * (len(self.levels))
+        for i in range(len(self.levels)):
+            print('level # {}, @ {}'.format(i, self.levels[i]))
+            lines = fn.isocurve(self.z, self.levels[i], connected=True, extendToEdge=True)[::-1]
+            lines = self._scale_contour_lines(lines)
+            oneline = self._join_lines(lines)
+            print('oneline ', oneline)
+
+            joined_lines[i] = oneline
+            x, y = map(list, zip(*oneline))
+            curves[i] = pg.PlotDataItem(x, y, pen=use_pen)#pens[i])
+            #if i == len(self.levels)-1:
+            #    self.ax.addItem(curve, pen=pg.mkPen(color='r'))
+            #if i == len(self.levels)-2:
+
+            # self.ax.addItem(curve)
+            # if i == 0:
+            #     x0 = np.mean([point[0] for point in oneline])
+            #     y0 = np.mean([point[1] for point in oneline])
+            #     xc = [x0, x0 + 1e-12]
+            #     yc = [y0, y0 + 1e-12]
+            #     curve_c = pg.PlotDataItem(xc, yc, pen=pens[i])
+            #     fill = pg.FillBetweenItem(curve, curve_c, brush=pg.mkBrush(color=self.colors[i]))
+            #     self.ax.addItem(fill)
+            # else:  # i > 0:
+            # #elif i == len(self.levels)-1:
+            #     fill = pg.FillBetweenItem(curve, prev_curve, brush=pg.mkBrush(color=self.colors[i]))
+            #     self.ax.addItem(fill)
+            # prev_curve = curve
+
+        # Get the curves at the edges of the array
+        x0 = np.mean([point[0] for point in joined_lines[1]])
+        y0 = np.mean([point[1] for point in joined_lines[1]])
+        dx0 = np.std([point[0] for point in joined_lines[1]])
+        dy0 = np.std([point[1] for point in joined_lines[1]])
+        x1 = np.mean([point[0] for point in joined_lines[-2]])
+        y1 = np.mean([point[1] for point in joined_lines[-2]])
+        dx1 = np.std([point[0] for point in joined_lines[-2]])
+        dy1 = np.std([point[1] for point in joined_lines[-2]])
+        if (dx1**2 + dy1**2) > (dx0**2 + dy0**2):
+            curves = [pg.PlotDataItem([x0, x0 + 1e-12], [y0, y0 + 1e-12], pen=pens[0])] + curves
+            #curves[-1] = pg.PlotDataItem(
+            #    np.array([self.x.min(), self.x.max()])[np.array([0, 0, 1, 1, 0])],
+            #    np.array([self.y.min(), self.y.max()])[np.array([0, 1, 1, 0, 0])], pen=pens[-1])
+        else:
+            curves = curves + [pg.PlotDataItem([x1, x1 + 1e-12], [y1, y1 + 1e-12], pen=pens[-1])]
+            #curves[0] = pg.PlotDataItem(
+            #    np.array([self.x.min(), self.x.max()])[0, 0, 1, 1, 0],
+            #    np.array([self.y.min(), self.y.max()])[0, 1, 1, 0, 0], pen=pens[0])
+
+        for j in range(len(self.levels)):
+            i = len(self.levels)-j-1
+            fill = pg.FillBetweenItem(curves[i], curves[i+1], brush=pg.mkBrush(color=self.colors[i]))
+            self.ax.addItem(fill)
+
     def draw_unfilled(self):
-        lws, lss = self.extl(self.linewidths), self.extl(self.linestyles)
+        lws, lss = self._extl(self.linewidths), self._extl(self.linestyles)
         pens = [setup_pen_kw(penkw=dict(color=self.colors[i]), linestyle=lss[i], linewidth=lws[i])
                 for i in range(len(self.levels))]
         contours = [pg.IsocurveItem(data=self.z, level=lvl, pen=pens[i]) for i, lvl in enumerate(self.levels)]
diff --git a/test.sh b/test.sh
index 195ad3d..310df8c 100755
--- a/test.sh
+++ b/test.sh
@@ -1,2 +1,2 @@
 #!/bin/bash
-python2.7 -m unittest discover --pattern=*.py -s tests
+python3.7 -m unittest discover --pattern=*.py -s tests
diff --git a/test3.6.sh b/test3.6.sh
new file mode 100755
index 0000000..0d855f5
--- /dev/null
+++ b/test3.6.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+# Script for launching tests in python 3 environment while otherwise working in python 2
+
+# Run the test
+export NEWPYTHONPATH=/lib/python/python3.6/site-packages/
+export PYTHONPATH=$NEWPYTHONPATH
+python3.6 -m unittest discover --pattern=*.py -s tests
+export PYTHONPATH=$OLDPYTHONPATH
+
+
+# More information
+# using pip with specific python versions:
+#     https://stackoverflow.com/a/33964956/6605826
+#     sudo python3.6 -m pip install -r requirements.txt
diff --git a/tests/examples.py b/tests/examples.py
index fd90f4e..e01bc5a 100755
--- a/tests/examples.py
+++ b/tests/examples.py
@@ -167,6 +167,20 @@ def twod_demo():
     return fig, axs
 
 
+def contour_demo():
+    """
+    Test contour-based plot methods
+    """
+    x = np.linspace(0, 1.8, 30)
+    y = np.linspace(1, 3.1, 25)
+    z = (x[:, np.newaxis] - 0.94)**2 + (y[np.newaxis, :] - 2.2)**2 + 1.145
+    levels = [1.2, 1.5, 2, 2.5, 2.95]
+    nlvl = len(levels) * 4
+
+    fig, axs = plt.subplots(2, 2)
+    return fig, axs
+
+
 def open_examples(close_after=False, start_event=True):
     print('pgmpl examples...')
     pgmpl.util.set_debug(0)
diff --git a/tests/test_axes.py b/tests/test_axes.py
index bf1b257..fb8367c 100755
--- a/tests/test_axes.py
+++ b/tests/test_axes.py
@@ -82,6 +82,9 @@ def test_axes_contour(self):
         levels = [0, 0.5, 1.2, 5, 9, 10, 20, 30]
         ax = Axes()
         ax.contour(a)
+        ax1 = Axes()
+        ax1.contourf(a)
+
         ax2 = Axes()
         ax2.contour(a, levels)
         ax3 = Axes()
diff --git a/tests/test_contour.py b/tests/test_contour.py
index a99f2a2..d34a0de 100755
--- a/tests/test_contour.py
+++ b/tests/test_contour.py
@@ -23,9 +23,21 @@ class TestPgmplContour(unittest.TestCase):
 
     verbose = int(os.environ.get('PGMPL_TEST_VERBOSE', '0'))
 
+    x0 = 0.94
+    y0 = 2.2
+    x02 = x0 - 0.45
+    y02 = y0 - 0.7
+    z0 = 1.145
+    slant = .0
     x = np.linspace(0, 1.8, 30)
     y = np.linspace(1, 3.1, 25)
-    z = (x[:, np.newaxis] - 0.94)**2 + (y[np.newaxis, :] - 2.2)**2 + 1.145
+    dx = x[np.newaxis, :] - x0
+    dy = y[:, np.newaxis] - y0
+    dx2 = x[np.newaxis, :] - x02
+    dy2 = y[:, np.newaxis] - y02
+    z = dx**2 + dy**2 + z0 + slant*(dx+dy)**2
+    z2 = -(dx2**2 + dy2**2 + (slant + 1.2) * (dx2 + dy2)**2)
+    z3 = dx**2 + dy**2 + (slant - 1.2) * (dx + dy)**2
     levels = [1.2, 1.5, 2, 2.5, 2.95]
     nlvl = len(levels) * 4
 
@@ -35,6 +47,7 @@ def printv(self, *args):
 
     def test_contour(self):
         fig, axs = subplots(4, 2)
+        fig.suptitle('TestPgmplContour.test_contour')
         axs[0, 0].set_title('z')
         axs[0, 0].contour(self.z)
         axs[0, 1].set_title('-z')
@@ -55,9 +68,96 @@ def test_contour(self):
         axs[3, 1].set_title('x, y, z, nlvl')
         axs[3, 1].contour(self.x, self.y, self.z, self.nlvl, linestyles=['-', '--', '-.', ':'])
 
+        # Repeat with matplotlib to check whether the same figure is produced
+        import matplotlib
+        from matplotlib import pyplot
+        fig, axs = pyplot.subplots(4, 2)
+        fig.suptitle('TestPgmplContour.test_contour')
+        axs[0, 0].set_title('z')
+        axs[0, 0].contour(self.z)
+        axs[0, 1].set_title('-z')
+        axs[0, 1].contour(-self.z, linestyles=['--', '-.', ':'])
+
+        axs[1, 0].set_title('z, levels')
+        axs[1, 0].contour(self.z, self.levels)
+        axs[1, 1].set_title('z, nlvl')
+        axs[1, 1].contour(self.z, self.nlvl, linewidths=[3, 2, 1])
+
+        axs[2, 0].set_title('x, y, z')
+        axs[2, 0].contour(self.x, self.y, self.z)
+        axs[2, 0].set_title('x, y, -z')
+        axs[2, 1].contour(self.x, self.y, -self.z, colors=['r', 'g', 'b'])
+
+        axs[3, 0].set_title('x, y, z, levels')
+        axs[3, 0].contour(self.x, self.y, self.z, self.levels)
+        axs[3, 1].set_title('x, y, z, nlvl')
+        axs[3, 1].contour(self.x, self.y, self.z, self.nlvl, linestyles=['-', '--', '-.', ':'])
+
+        # # Uncomment these to see figures during manual testing and development:
+        # pyplot.show()
+        # import pgmpl
+        # pgmpl.app.exec_()
+
     def test_contourf(self):
-        fig, ax = subplots(1)
-        self.assertRaises(NotImplementedError, ax.contourf, self.z)
+        fig, axs = subplots(4, 2)
+        fig.suptitle('TestPgmplContour.test_contourf')
+        axs[0, 0].set_title('z')
+        axs[0, 0].contourf(self.z)
+
+        axs[0, 1].set_title('z, nlvl')
+        axs[0, 1].contourf(self.z, self.nlvl)
+
+        axs[1, 0].set_title('z2')
+        axs[1, 0].contourf(self.z2)
+
+        axs[1, 1].set_title('z3')
+        axs[1, 1].contourf(self.z3)
+
+        axs[2, 0].set_title('x, y, z')
+        axs[2, 0].contourf(self.x, self.y, self.z)
+
+        axs[2, 1].set_title('x, y, z, nlvl')
+        axs[2, 1].contourf(self.x, self.y, self.z, self.nlvl)
+
+        axs[3, 0].set_title('x, y, z2')
+        axs[3, 0].contourf(self.x, self.y, self.z2)
+
+        axs[3, 1].set_title('x, y, z3')
+        axs[3, 1].contourf(self.x, self.y, self.z3)
+
+        # Repeat with matplotlib to check whether the same figure is produced
+        import matplotlib
+        from matplotlib import pyplot
+        figm, axsm = pyplot.subplots(4, 2)
+        figm.suptitle('TestPgmplContour.test_contourf')
+        axsm[0, 0].set_title('z')
+        axsm[0, 0].contourf(self.z)
+
+        axsm[0, 1].set_title('z, nlvl')
+        axsm[0, 1].contourf(self.z, self.nlvl)
+
+        axsm[1, 0].set_title('z2')
+        axsm[1, 0].contourf(self.z2)
+
+        axsm[1, 1].set_title('z3')
+        axsm[1, 1].contourf(self.z3)
+
+        axsm[2, 0].set_title('x, y, z')
+        axsm[2, 0].contourf(self.x, self.y, self.z)
+
+        axsm[2, 1].set_title('x, y, z, nlvl')
+        axsm[2, 1].contourf(self.x, self.y, self.z, self.nlvl)
+
+        axsm[3, 0].set_title('x, y, z2')
+        axsm[3, 0].contourf(self.x, self.y, self.z2)
+
+        axsm[3, 1].set_title('x, y, z3')
+        axsm[3, 1].contourf(self.x, self.y, self.z3)
+
+        # Uncomment these to see figures during manual testing and development:
+        pyplot.show()
+        import pgmpl
+        pgmpl.app.exec_()
 
     def test_contour_errors(self):
         ax = Axes()