diff --git a/.pylintrc b/.pylintrc
index 66b76a52b7c..9e16ab9da39 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -441,6 +441,7 @@ function-naming-style=snake_case
 good-names=i,
            j,
            k,
+           ax,
            ex,
            Run,
            _,
diff --git a/doc/api/index.rst b/doc/api/index.rst
index c4f199b3dbc..3ee226b9afc 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -16,6 +16,7 @@ All plotting is handled through the :class:`pygmt.Figure` class and its methods.
     :toctree: generated
 
     Figure
+    subplots
 
 Plotting data and laying out the map:
 
diff --git a/doc/index.rst b/doc/index.rst
index 4c3fd3d34e3..33a005631f2 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -32,6 +32,7 @@
     projections/index.rst
     tutorials/coastlines.rst
     tutorials/plot.rst
+    tutorials/subplots.rst
     tutorials/text.rst
 
 .. toctree::
diff --git a/examples/tutorials/subplots.py b/examples/tutorials/subplots.py
new file mode 100644
index 00000000000..a6fcef21abf
--- /dev/null
+++ b/examples/tutorials/subplots.py
@@ -0,0 +1,224 @@
+"""
+Subplots
+========
+
+When you're preparing a figure for a paper, there will often be times when
+you'll need to put many individual plots into one large figure, and label them
+'abcd'. These individual plots are called subplots.
+
+There are two main ways to create subplots in GMT:
+
+- Use :meth:`pygmt.Figure.shift_origin` to manually move each individual plot
+  to the right position.
+- Use :meth:`pygmt.subplots` to define the layout of the subplots.
+
+The first method is easier to use and should handle simple cases involving a
+couple of subplots. For more advanced subplot layouts however, we recommend the
+use of :meth:`pygmt.subplots` which offers finer grained control, and this is
+what the tutorial below will cover.
+"""
+
+###############################################################################
+# Let's start by importing the PyGMT library
+
+import pygmt
+
+###############################################################################
+# Define subplot layout
+# ---------------------
+#
+# The :meth:`pygmt.subplots` command is used to setup the layout, size, and
+# other attributes of the figure. It divides the whole canvas into regular grid
+# areas with n rows and m columns. Each grid area can contain an individual
+# subplot. For example:
+
+fig, axs = pygmt.subplots(nrows=2, ncols=3, figsize=("15c", "6c"), frame="lrtb")
+
+###############################################################################
+# will define our figure to have a 2 row and 3 column grid layout.
+# ``figsize=("15c", "6c")`` defines the overall size of the figure to be 15cm
+# wide by 6cm high. Using ``frame="lrtb"`` allows us to customize the map frame
+# for all subplots instead of setting them individually. The figure layout will
+# look like the following:
+
+for index in axs.flatten():
+    i = index // axs.shape[1]  # row
+    j = index % axs.shape[1]  # column
+    fig.sca(ax=axs[i, j])  # sets the current Axes
+    fig.text(
+        position="MC", text=f"index: {index}, row: {i}, col: {j}", region=[0, 1, 0, 1]
+    )
+fig.end_subplot()
+fig.show()
+
+###############################################################################
+# The ``fig.sca`` command activates a specified subplot, and all subsequent
+# plotting commands will take place in that subplot. This is similar to
+# matplotlib's ``plt.sca`` method. In order to specify a subplot, you will need
+# to provide the identifier for that subplot via the ``ax`` argument. This can
+# be found in the ``axs`` variable referenced by the ``row`` and ``col``
+# number.
+
+###############################################################################
+# .. note::
+#
+#     The row and column numbering starts from 0. So for a subplot layout with
+#     N rows and M columns, row numbers will go from 0 to N-1, and column
+#     numbers will go from 0 to M-1.
+
+###############################################################################
+# For example, to activate the subplot on the top right corner (index: 2) at
+# ``row=0`` and ``col=2``, so that all subsequent plotting commands happen
+# there, you can use the following command:
+
+###############################################################################
+# .. code-block:: default
+#
+#     fig.sca(ax=axs[0, 2])
+
+###############################################################################
+# Finally, remember to use ``fig.end_subplot()`` to exit the subplot mode.
+
+###############################################################################
+# .. code-block:: default
+#
+#     fig.end_subplot()
+
+###############################################################################
+# Making your first subplot
+# -------------------------
+# Next, let's use what we learned above to make a 2 row by 2 column subplot
+# figure. We'll also pick up on some new parameters to configure our subplot.
+
+fig, axs = pygmt.subplots(
+    nrows=2,
+    ncols=2,
+    figsize=("15c", "6c"),
+    autolabel=True,
+    margins=["0.1c", "0.2c"],
+    title="My Subplot Heading",
+)
+fig.basemap(region=[0, 10, 0, 10], projection="X?", frame=["af", "WSne"], ax=axs[0, 0])
+fig.basemap(region=[0, 20, 0, 10], projection="X?", frame=["af", "WSne"], ax=axs[0, 1])
+fig.basemap(region=[0, 10, 0, 20], projection="X?", frame=["af", "WSne"], ax=axs[1, 0])
+fig.basemap(region=[0, 20, 0, 20], projection="X?", frame=["af", "WSne"], ax=axs[1, 1])
+fig.end_subplot()
+fig.show()
+
+###############################################################################
+# In this example, we define a 2-row, 2-column (2x2) subplot layout using
+# :meth:`pygmt.subplots`. The overall figure dimensions is set to be 15cm wide
+# and 6cm high (``figsize=["15c", "6c"]``). In addition, we used some optional
+# parameters to fine tune some details of the figure creation:
+#
+# - ``autolabel=True``: Each subplot is automatically labelled abcd
+# - ``margins=["0.1c", "0.2c"]``: adjusts the space between adjacent subplots.
+#   In this case, it is set as 0.1 cm in the X direction and 0.2 cm in the Y
+#   direction.
+# - ``title="My Subplot Heading"``: adds a title on top of the whole figure.
+#
+# Notice that each subplot was set to use a linear projection ``"X?"``.
+# Usually, we need to specify the width and height of the map frame, but it is
+# also possible to use a question mark ``"?"`` to let GMT decide automatically
+# on what is the most appropriate width/height for the each subplot's map
+# frame.
+
+###############################################################################
+# .. tip::
+#
+#     In the above example, we used the following commands to activate the
+#     four subplots explicitly one after another::
+#
+#         fig.basemap(..., ax=axs[0, 0])
+#         fig.basemap(..., ax=axs[0, 1])
+#         fig.basemap(..., ax=axs[1, 0])
+#         fig.basemap(..., ax=axs[1, 1])
+#
+#     In fact, we can just use ``fig.basemap(..., ax=True)`` without specifying
+#     any subplot index number, and GMT will automatically activate the next
+#     subplot.
+
+###############################################################################
+# Shared X and Y axis labels
+# --------------------------
+# In the example above with the four subplots, the two subplots for each row
+# have the same Y-axis range, and the two subplots for each column have the
+# same X-axis range. You can use the **layout** option to set a common X and/or
+# Y axis between subplots.
+
+fig, axs = pygmt.subplots(
+    nrows=2,
+    ncols=2,
+    figsize=("15c", "6c"),
+    autolabel=True,
+    margins=["0.3c", "0.2c"],
+    title="My Subplot Heading",
+    layout=["Rl", "Cb"],
+    frame="WSrt",
+)
+fig.basemap(region=[0, 10, 0, 10], projection="X?", ax=True)
+fig.basemap(region=[0, 20, 0, 10], projection="X?", ax=True)
+fig.basemap(region=[0, 10, 0, 20], projection="X?", ax=True)
+fig.basemap(region=[0, 20, 0, 20], projection="X?", ax=True)
+fig.end_subplot()
+fig.show()
+
+###############################################################################
+# **Rl** indicates that subplots within a **R**\ ow will share the y-axis, and
+# only the **l**\ eft axis is displayed. **Cb** indicates that subplots in
+# a column will share the x-axis, and only the **b**\ ottom axis is displayed.
+#
+# Of course, instead of using the **layout** option, you can also set a
+# different **frame** for each subplot to control the axis properties
+# individually for each subplot.
+
+###############################################################################
+# Advanced subplot layouts
+# ------------------------
+#
+# Nested subplot are currently not supported. If you want to create more
+# complex subplot layouts, some manual adjustments are needed.
+#
+# The following example draws three subplots in a 2-row, 2-column layout, with
+# the first subplot occupying the first row.
+
+fig, axs = pygmt.subplots(nrows=2, ncols=2, figsize=("15c", "6c"), autolabel=True)
+fig.basemap(
+    region=[0, 10, 0, 10], projection="X15c/3c", frame=["af", "WSne"], ax=axs[0, 0]
+)
+fig.text(text="TEXT", x=5, y=5, projection="X15c/3c")
+fig.basemap(region=[0, 5, 0, 5], projection="X?", frame=["af", "WSne"], ax=axs[1, 0])
+fig.basemap(region=[0, 5, 0, 5], projection="X?", frame=["af", "WSne"], ax=axs[1, 1])
+fig.end_subplot()
+fig.show()
+
+###############################################################################
+#
+# When drawing the three basemaps, the last two basemaps use
+# ``projection="X?"``, so GMT will automatically determine the size of the
+# subplot according to the size of the subplot area. In order for the first
+# subplot to fill up the entire top row space, we use manually adjusted the
+# subplot width to 15cm using ``projection="X15c/3c"``.
+
+###############################################################################
+# .. note::
+#
+#     There are bugs that have not been fixed in the above example.
+#
+#     In subplot mode, the size of each subgraph is controlled by the
+#     ``figsize`` option of :meth:`pygmt.subplots`. Users can override this and
+#     use``projection`` to specify the size of an individual subplot, but this
+#     size will not be remembered. If the next command does not specify
+#     ``projection``, the default size of the subplot mode will be used, and
+#     the resulting plot will be inccorect.
+#
+#     The current workaround is to use the same ``projection`` option in all
+#     commands for the subplot. For example, we forced subplot (a) to have a
+#     different size using ``projection="15c/3c``. The next command within the
+#     subplot (e.g. ``text``) must also use ``projection="x15c/3c"``, otherwise
+#     the placement will be wrong.
+
+###############################################################################
+# Since we skipped the second subplot, the auto label function will name the
+# three subplots as a, c and d, which is not what we want, so we have to use
+# ``fig.sca(A=""(a)"`` to manually set the subplot label.
diff --git a/pygmt/__init__.py b/pygmt/__init__.py
index 512865264b0..ef7792824d2 100644
--- a/pygmt/__init__.py
+++ b/pygmt/__init__.py
@@ -19,6 +19,7 @@
 from .sampling import grdtrack
 from .mathops import makecpt
 from .modules import GMTDataArrayAccessor, config, info, grdinfo, which
+from .subplot import SubPlot, subplots
 from .gridops import grdcut
 from .x2sys import x2sys_init, x2sys_cross
 from . import datasets
diff --git a/pygmt/base_plotting.py b/pygmt/base_plotting.py
index 0f7a20ad7f3..3f1720c9209 100644
--- a/pygmt/base_plotting.py
+++ b/pygmt/base_plotting.py
@@ -70,6 +70,7 @@ def _preprocess(self, **kwargs):  # pylint: disable=no-self-use
         U="timestamp",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
@@ -153,6 +154,7 @@ def coast(self, **kwargs):
         W="scale",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
@@ -238,9 +240,10 @@ def colorbar(self, **kwargs):
         S="resample",
         U="timestamp",
         W="pen",
-        l="label",
         X="xshift",
         Y="yshift",
+        c="ax",
+        l="label",
         p="perspective",
         t="transparency",
     )
@@ -334,6 +337,7 @@ def grdcontour(self, grid, **kwargs):
         V="verbose",
         X="xshift",
         Y="yshift",
+        c="ax",
         n="interpolation",
         p="perspective",
         t="transparency",
@@ -490,6 +494,7 @@ def grdimage(self, grid, **kwargs):
         I="shading",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
@@ -609,6 +614,7 @@ def grdview(self, grid, **kwargs):
         U="timestamp",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
@@ -735,6 +741,7 @@ def plot(self, x=None, y=None, data=None, sizes=None, direction=None, **kwargs):
         C="levels",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
@@ -829,6 +836,7 @@ def contour(self, x=None, y=None, z=None, data=None, **kwargs):
         U="timestamp",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
@@ -870,7 +878,7 @@ def basemap(self, **kwargs):
 
         """
         kwargs = self._preprocess(**kwargs)
-        if not ("B" in kwargs or "L" in kwargs or "T" in kwargs):
+        if not ("B" in kwargs or "L" in kwargs or "T" in kwargs or "c" in kwargs):
             raise GMTInvalidInput("At least one of B, L, or T must be specified.")
         with Session() as lib:
             lib.call_module("basemap", build_arg_string(kwargs))
@@ -884,6 +892,7 @@ def basemap(self, **kwargs):
         F="box",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
@@ -932,6 +941,7 @@ def logo(self, **kwargs):
         M="monochrome",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
@@ -985,6 +995,7 @@ def image(self, imagefile, **kwargs):
         F="box",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
@@ -1055,6 +1066,7 @@ def legend(self, spec=None, position="JTR+jTR+o0.2c", box="+gwhite+p1p", **kwarg
         W="pen",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
@@ -1222,6 +1234,7 @@ def text(
         C="offset",
         X="xshift",
         Y="yshift",
+        c="ax",
         p="perspective",
         t="transparency",
     )
diff --git a/pygmt/subplot.py b/pygmt/subplot.py
new file mode 100644
index 00000000000..7c92232ea5c
--- /dev/null
+++ b/pygmt/subplot.py
@@ -0,0 +1,228 @@
+"""High level functions for making subplots."""
+import numpy as np
+
+from .clib import Session
+from .figure import Figure
+from .helpers import build_arg_string, fmt_docstring, kwargs_to_strings, use_alias
+
+
+class SubPlot(Figure):
+    """
+    Manage modern mode figure subplot configuration and selection.
+
+    The subplot module is used to split the current figure into a
+    rectangular layout of subplots that each may contain a single
+    self-contained figure. A subplot setup is started with the begin
+    directive that defines the layout of the subplots, while positioning to
+    a particular subplot for plotting is done via the set directive. The
+    subplot process is completed via the end directive.
+
+    Full option list at :gmt-docs:`subplot.html`
+    """
+
+    def __init__(self, nrows, ncols, figsize, **kwargs):
+        super().__init__()
+        # Activate main Figure, and initiate subplot
+        self._activate_figure()
+        self.begin_subplot(row=nrows, col=ncols, figsize=figsize, **kwargs)
+
+    @staticmethod
+    @fmt_docstring
+    @use_alias(
+        Ff="figsize",
+        A="autolabel",
+        B="frame",
+        C="clearance",
+        M="margins",
+        S="layout",
+        T="title",
+    )
+    @kwargs_to_strings(Ff="sequence", M="sequence")
+    def begin_subplot(row=None, col=None, **kwargs):
+        """
+        The begin directive of subplot defines the layout of the entire multi-
+        panel illustration. Several options are available to specify the
+        systematic layout, labeling, dimensions, and more for the subplots.
+
+        {aliases}
+        """
+        arg_str = " ".join(["begin", f"{row}x{col}", build_arg_string(kwargs)])
+        with Session() as lib:
+            lib.call_module(module="subplot", args=arg_str)
+
+    @staticmethod
+    @fmt_docstring
+    @use_alias(F="dimensions")
+    def sca(ax=None, **kwargs):
+        """
+        Set the current Axes instance to *ax*.
+
+        Before you start plotting you must first select the active subplot.
+        Note: If any projection (J) option is passed with ? as scale or
+        width when plotting subplots, then the dimensions of the map are
+        automatically determined by the subplot size and your region. For
+        Cartesian plots: If you want the scale to apply equally to both
+        dimensions then you must specify ``projection="x"`` [The default
+        ``projection="X"`` will fill the subplot by using unequal scales].
+
+        {aliases}
+        """
+        arg_str = " ".join(["set", f"{ax}", build_arg_string(kwargs)])
+        with Session() as lib:
+            lib.call_module(module="subplot", args=arg_str)
+
+    @staticmethod
+    @fmt_docstring
+    @use_alias(V="verbose")
+    def end_subplot(**kwargs):
+        """
+        This command finalizes the current subplot, including any placement of
+        tags, and updates the gmt.history to reflect the dimensions and linear
+        projection required to draw the entire figure outline. This allows
+        subsequent commands, such as colorbar, to use ``position="J"`` to place
+        bars with reference to the complete figure dimensions. We also reset
+        the current plot location to where it was prior to the subplot.
+
+        {aliases}
+        """
+        arg_str = " ".join(["end", build_arg_string(kwargs)])
+        with Session() as lib:
+            lib.call_module(module="subplot", args=arg_str)
+
+
+def subplots(
+    nrows=1,
+    ncols=1,
+    figsize=(6.4, 4.8),
+    autolabel=None,
+    clearance=None,
+    margins=None,
+    layout=None,
+    title=None,
+    **kwargs,
+):
+    """
+    Create a figure with a set of subplots.
+
+    Full option list at :gmt-docs:`subplot.html#synopsis-begin-mode`
+
+    Parameters
+    ----------
+    nrows : int
+        Number of rows of the subplot grid.
+
+    ncols : int
+        Number of columns of the subplot grid.
+
+    figsize : tuple
+        Overall figure dimensions as ``(width, height)``. Default is (6.4, 4.8)
+
+    autolabel : bool or str
+        ``[autolabel][+cdx[/dy]][+gfill][+j|Jrefpoint][+odx[/dy]][+ppen][+r|R]
+        [+v]``.
+        Specify automatic tagging of each subplot. Append either a number or
+        letter [a]. This sets the tag of the first, top-left subplot and others
+        follow sequentially. Surround the number or letter by parentheses on
+        any side if these should be typeset as part of the tag. Use
+        **+j|J**\\ *refpoint* to specify where the tag should be placed in the
+        subplot [TL]. Note: **+j** sets the justification of the tag to
+        *refpoint* (suitable for interior tags) while **+J** instead selects
+        the mirror opposite (suitable for exterior tags). Append
+        **+c**\\ *dx*[/*dy*] to set the clearance between the tag and a
+        surrounding text box requested via **+g** or **+p** [3p/3p, i.e., 15%
+        of the FONT_TAG size dimension]. Append **+g**\\ *fill* to paint the
+        tag's text box with *fill* [no painting]. Append
+        **+o**\\ *dx*\\ [/*dy*] to offset the tag's reference point in the
+        direction implied by the justification [4p/4p, i.e., 20% of the
+        FONT_TAG size]. Append **+p**\\ *pen* to draw the outline of the tag's
+        text box using selected *pen* [no outline]. Append **+r** to typeset
+        your tag numbers using lowercase Roman numerals; use **+R** for
+        uppercase Roman numerals [Arabic numerals]. Append **+v** to increase
+        tag numbers vertically down columns [horizontally across rows].
+
+    clearance : str
+        ``[side]clearance``.
+        Reserve a space of dimension *clearance* between the margin and the
+        subplot on the specified side, using *side* values from **w**, **e**,
+        **s**, or **n**, or **x** for both **w** and **e** or **y** for both
+        **s** and **n**.  No *side* means all sides. The option is repeatable
+        to set aside space on more than one side. Such space will be left
+        untouched by the main map plotting but can be accessed by modules that
+        plot scales, bars, text, etc.  Settings specified under **begin**
+        directive apply to all subplots, while settings under **set** only
+        apply to the selected (active) subplot. **Note**: Common options
+        **x_offset** and **y_offset* are not available during subplots; use
+        **clearance** instead.
+
+    margins : tuple
+        This is margin space that is added between neighboring subplots (i.e.,
+        the interior margins) in addition to the automatic space added for tick
+        marks, annotations, and labels. The margins can be specified as either:
+
+        - a single value (for same margin on all sides). E.g. '5c'.
+        - a pair of values (for setting separate horizontal and vertical
+          margins). E.g. ['5c', '3c'].
+        - a set of four values (for setting separate left, right, bottom, and
+          top margins). E.g. ['1c', '2c', '3c', '4c'].
+
+        The actual gap created is always a sum of the margins for the two
+        opposing sides (e.g., east plus west or south plus north margins)
+        [Default is half the primary annotation font size, giving the full
+        annotation font size as the default gap].
+
+    layout : str or list
+        Set subplot layout for shared axes. May be set separately for rows
+        (**R**) and columns (**C**). E.g. ``layout=['Rl', 'Cb']``.
+        Considerations for **C**: Use when all subplots in a **C**\\ olumn
+        share a common *x*-range. The first (i.e., **t**\\ op) and the last
+        (i.e., **b**\\ ottom) rows will have *x* annotations; append **t** or
+        **b** to select only one of those two rows [both]. Append **+l** if
+        annotated *x*-axes should have a label [none]; optionally append the
+        label if it is the same for the entire subplot. Append **+t** to make
+        space for subplot titles for each row; use **+tc** for top row titles
+        only [no subplot titles]. Labels and titles that depends on which row
+        or column are specified as usual via a subplot's own **frame** setting.
+        Considerations for **R**: Use when all subplots in a **R**\\ ow share a
+        common *y*-range. The first (i.e., **l**\\ eft) and the last (i.e.,
+        **r**\\ ight) columns will have *y*-annotations; append **l** or **r**
+        to select only one of those two columns [both]. Append **+l** if
+        annotated *y*-axes will have a label [none]; optionally, append the
+        label if it is the same for the entire subplot. Append **+p** to make
+        all annotations axis-parallel [horizontal]; if not used you may have to
+        set **clearance** to secure extra space for long horizontal
+        annotations. Append **+w** to draw horizontal and vertical lines
+        between interior panels using selected pen [no lines].
+
+    title : str
+        Overarching heading for the entire figure. Font is determined by
+        setting ``FONT_HEADING``.
+
+    Returns
+    -------
+    fig : :class:`pygmt.Figure`
+        A PyGMT Figure instance.
+
+    axs : numpy.ndarray
+        Array of Axes objects.
+    """
+    # Get PyGMT Figure with SubPlot initiated
+    fig = SubPlot(
+        nrows=nrows,
+        ncols=ncols,
+        figsize=figsize,
+        autolabel=autolabel,
+        clearance=clearance,
+        margins=margins,
+        layout=layout,
+        title=f'"{title}"' if title else None,
+        **kwargs,
+    )
+
+    # Setup matplotlib-like Axes
+    axs = np.empty(shape=(nrows, ncols), dtype=object)
+    for index in range(nrows * ncols):
+        i = index // ncols  # row
+        j = index % ncols  # column
+        axs[i, j] = index
+
+    return fig, axs
diff --git a/pygmt/tests/baseline/test_subplot_basic.png b/pygmt/tests/baseline/test_subplot_basic.png
new file mode 100644
index 00000000000..6a9a6a6b5e5
Binary files /dev/null and b/pygmt/tests/baseline/test_subplot_basic.png differ
diff --git a/pygmt/tests/baseline/test_subplot_direct.png b/pygmt/tests/baseline/test_subplot_direct.png
new file mode 100644
index 00000000000..4cedf3e8669
Binary files /dev/null and b/pygmt/tests/baseline/test_subplot_direct.png differ
diff --git a/pygmt/tests/baseline/test_subplot_frame.png b/pygmt/tests/baseline/test_subplot_frame.png
new file mode 100644
index 00000000000..071133e1869
Binary files /dev/null and b/pygmt/tests/baseline/test_subplot_frame.png differ
diff --git a/pygmt/tests/test_subplot.py b/pygmt/tests/test_subplot.py
new file mode 100644
index 00000000000..50ee378d99f
--- /dev/null
+++ b/pygmt/tests/test_subplot.py
@@ -0,0 +1,90 @@
+"""
+Tests subplot
+"""
+import pytest
+
+from ..helpers.testing import check_figures_equal
+from ..subplot import SubPlot, subplots
+
+
+@pytest.mark.mpl_image_compare
+def test_subplot_basic():
+    """
+    Create a subplot figure with 1 row and 2 columns.
+    """
+    fig, axs = subplots(nrows=1, ncols=2, figsize=("6c", "3c"))
+    fig.sca(ax=axs[0, 0])
+    fig.basemap(region=[0, 3, 0, 3], frame=True)
+    fig.sca(ax=axs[0, 1])
+    fig.basemap(region=[0, 3, 0, 3], frame=True)
+    fig.end_subplot()
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+def test_subplot_frame():
+    """
+    Check that map frame setting is applied to all subplot figures
+    """
+    fig, axs = subplots(nrows=1, ncols=2, figsize=("6c", "3c"), frame="WSne")
+    fig.sca(ax=axs[0, 0])
+    fig.basemap(region=[0, 3, 0, 3], frame="+tplot0")
+    fig.sca(ax=axs[0, 1])
+    fig.basemap(region=[0, 3, 0, 3], frame="+tplot1")
+    fig.end_subplot()
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+def test_subplot_direct():
+    """
+    Plot map elements to subplots directly using ax argument
+    """
+    fig, axs = subplots(nrows=2, ncols=1, figsize=("3c", "6c"))
+    fig.basemap(region=[0, 3, 0, 3], frame=True, ax=axs[0, 0])
+    fig.basemap(region=[0, 3, 0, 3], frame=True, ax=axs[1, 0])
+    fig.end_subplot()
+    return fig
+
+
+@check_figures_equal()
+def test_subplot_autolabel_margins_title():
+    """
+    Make subplot figure with autolabels, setting some margins and a title.
+    """
+    kwargs = dict(nrows=2, ncols=1, figsize=("15c", "6c"))
+
+    fig_ref = SubPlot(A="(1)", M="0.3c/0.1c", T='"Subplot Title"', **kwargs)
+    fig_ref.basemap(region=[0, 1, 2, 3], frame="WSne", c="0,0")
+    fig_ref.basemap(region=[4, 5, 6, 7], frame="WSne", c="1,0")
+    fig_ref.end_subplot()
+
+    fig_test, axs_test = subplots(
+        autolabel="(1)", margins=["0.3c", "0.1c"], title="Subplot Title", **kwargs
+    )
+    fig_test.basemap(region=[0, 1, 2, 3], frame="WSne", ax=axs_test[0, 0])
+    fig_test.basemap(region=[4, 5, 6, 7], frame="WSne", ax=axs_test[1, 0])
+    fig_test.end_subplot()
+
+    return fig_ref, fig_test
+
+
+@check_figures_equal()
+def test_subplot_clearance_and_shared_xy_axis_layout():
+    """
+    Ensure subplot clearance works, and that the layout can be set to use
+    shared X and Y axis labels across columns and rows.
+    """
+    kwargs = dict(nrows=2, ncols=2, frame="WSrt", figsize=("5c", "5c"))
+
+    fig_ref = SubPlot(C="y0.2", SR="l", SC="t", **kwargs)
+    fig_test, _ = subplots(clearance="y0.2", layout=["Rl", "Ct"], **kwargs)
+
+    for fig in (fig_ref, fig_test):
+        fig.basemap(region=[0, 4, 0, 4], projection="X?", ax=True)
+        fig.basemap(region=[0, 8, 0, 4], projection="X?", ax=True)
+        fig.basemap(region=[0, 4, 0, 8], projection="X?", ax=True)
+        fig.basemap(region=[0, 8, 0, 8], projection="X?", ax=True)
+        fig.end_subplot()
+
+    return fig_ref, fig_test