diff --git a/doc/api/index.rst b/doc/api/index.rst
index 0f673b90d5b..3b4220b7e7d 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -36,6 +36,7 @@ Plotting data and laying out the map:
     Figure.meca
     Figure.plot
     Figure.plot3d
+    Figure.rose
     Figure.set_panel
     Figure.shift_origin
     Figure.solar
diff --git a/examples/gallery/embellishments/rose.py b/examples/gallery/embellishments/rose.py
new file mode 100644
index 00000000000..1c4ca687ec6
--- /dev/null
+++ b/examples/gallery/embellishments/rose.py
@@ -0,0 +1,44 @@
+"""
+Rose diagram
+------------
+
+The :meth:`pygmt.Figure.rose` method can plot windrose diagrams or polar histograms.
+"""
+
+import pygmt
+
+# Load sample compilation of fracture lengths and azimuth as
+# hypothetically digitized from geological maps
+data = pygmt.datasets.load_fractures_compilation()
+
+fig = pygmt.Figure()
+
+fig.rose(
+    # use columns of the sample dataset as input for the length and azimuth
+    # parameters
+    length=data.length,
+    azimuth=data.azimuth,
+    # specify the "region" of interest in the (r,azimuth) space [r0, r1, az0, az1],
+    # here, r0 is 0 and r1 is 1, for azimuth, az0 is 0 and az1 is 360 which means
+    # we plot a full circle between 0 and 360 degrees
+    region=[0, 1, 0, 360],
+    # set the diameter of the rose diagram to 7.5 cm
+    diameter="7.5c",
+    # define the sector width in degrees, we append +r here to draw a rose
+    # diagram instead of a sector diagram
+    sector="10+r",
+    # normalize bin counts by the largest value so all bin counts range from
+    # 0 to 1
+    norm=True,
+    # use red3 as color fill for the sectors
+    color="red3",
+    # define the frame with ticks and gridlines every 0.2
+    # length unit in radial direction and every 30 degrees
+    # in azimuthal direction, set background color to
+    # lightgray
+    frame=["x0.2g0.2", "y30g30", "+glightgray"],
+    # use a pen size of 1p to draw the outlines
+    pen="1p",
+)
+
+fig.show()
diff --git a/pygmt/figure.py b/pygmt/figure.py
index a8d89b1e727..db152e0d27c 100644
--- a/pygmt/figure.py
+++ b/pygmt/figure.py
@@ -390,6 +390,7 @@ def _repr_html_(self):
         meca,
         plot,
         plot3d,
+        rose,
         set_panel,
         solar,
         subplot,
diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py
index 7fae8fd8a2d..7911195bf57 100644
--- a/pygmt/src/__init__.py
+++ b/pygmt/src/__init__.py
@@ -2,6 +2,7 @@
 Source code for PyGMT modules.
 """
 # pylint: disable=import-outside-toplevel
+
 from pygmt.src.basemap import basemap
 from pygmt.src.blockm import blockmean, blockmedian
 from pygmt.src.coast import coast
@@ -25,6 +26,7 @@
 from pygmt.src.meca import meca
 from pygmt.src.plot import plot
 from pygmt.src.plot3d import plot3d
+from pygmt.src.rose import rose
 from pygmt.src.solar import solar
 from pygmt.src.subplot import set_panel, subplot
 from pygmt.src.surface import surface
diff --git a/pygmt/src/rose.py b/pygmt/src/rose.py
new file mode 100644
index 00000000000..cd91b3d0799
--- /dev/null
+++ b/pygmt/src/rose.py
@@ -0,0 +1,204 @@
+"""
+rose - Plot windrose diagrams or polar histograms.
+"""
+
+from pygmt.clib import Session
+from pygmt.helpers import build_arg_string, fmt_docstring, kwargs_to_strings, use_alias
+
+
+@fmt_docstring
+@use_alias(
+    A="sector",
+    B="frame",
+    C="cmap",
+    D="shift",
+    Em="vectors",
+    F="no_scale",
+    G="color",
+    I="inquire",
+    JX="diameter",
+    L="labels",
+    M="vector_params",
+    Q="alpha",
+    R="region",
+    S="norm",
+    T="orientation",
+    U="timestamp",
+    V="verbose",
+    W="pen",
+    X="xshift",
+    Y="yshift",
+    Z="scale",
+    i="columns",
+    c="panel",
+    p="perspective",
+    t="transparency",
+)
+@kwargs_to_strings(R="sequence", c="sequence_comma", i="sequence_comma", p="sequence")
+def rose(self, length=None, azimuth=None, data=None, **kwargs):
+    """
+    Plot windrose diagrams or polar histograms.
+
+    Takes a matrix, (length,azimuth) pairs, or a file name as input
+    and plots windrose diagrams or polar histograms (sector diagram
+    or rose diagram).
+
+    Must provide either ``data`` or ``length`` and ``azimuth``.
+
+    Options include full circle and half circle plots. The outline
+    of the windrose is drawn with the same color as
+    :gmt-term:`MAP_DEFAULT_PEN`.
+
+    Full option list at :gmt-docs:`rose.html`
+
+    {aliases}
+
+    Parameters
+    ----------
+    length/azimuth : float or 1d arrays
+        Length and azimuth values, or arrays of length and azimuth
+        values
+
+    data : str or 2d array
+        Either a data file name or a 2d numpy array with the tabular data.
+        Use option ``columns`` to choose which columns are length and
+        azimuth, respectively. If a file with only azimuths is given,
+        use ``columns`` to indicate the single column with azimuths; then
+        all lengths are set to unity (see ``scale = 'u'`` to set actual
+        lengths to unity as well).
+
+    orientation : bool
+        Specifies that the input data are orientation data (i.e., have a
+        180 degree ambiguity) instead of true 0-360 degree directions
+        [Default is 0-360 degrees]. We compensate by counting each record
+        twice: First as azimuth and second as azimuth +180. Ignored if
+        range is given as -90/90 or 0/180.
+
+    region : str or list
+        *Required if this is the first plot command*.
+        *r0/r1/az0/az1*.
+        Specifies the 'region' of interest in (*r*, *azimuth*) space. Here,
+        *r0* is 0, *r1* is max length in units. For *az0* and *az1*,
+        specify either -90/90 or 0/180 for half circle plot or 0/360 for
+        full circle.
+
+    diameter : str
+         Sets the diameter of the rose diagram. If not given,
+         then we default to a diameter of 7.5 cm.
+
+    sector : str
+         Gives the sector width in degrees for sector and rose diagram.
+         Default ``'0'`` means windrose diagram. Append **+r** to draw rose
+         diagram instead of sector diagram (e.g. ``'10+r'``).
+
+    norm : bool
+         Normalize input radii (or bin counts if ``sector_width`` is used)
+         by the largest value so all radii (or bin counts) range from 0
+         to 1.
+
+    frame : str
+         Set map boundary frame and axes attributes. Remember that *x*
+         here is radial distance and *y* is azimuth. The ylabel may be
+         used to plot a figure caption. The scale bar length is determined
+         by the radial gridline spacing.
+
+    scale : float or str
+         Multiply the data radii by scale. E.g., use ``scale = 0.001`` to
+         convert your data from m to km. To exclude the radii from
+         consideration, set them all to unity with ``scale = 'u'``
+         [Default is no scaling].
+
+    columns : str or 1d array
+         Select input columns and transformations. E.g. choose
+         ``columns = [1, 0]`` or ``columns = '1,0'`` if the length values
+         are stored in the second column and the direction (azimuth)
+         values in the first one. Note: zero-based indexing is used.
+
+    color : str
+         Selects shade, color or pattern for filling the sectors [Default
+         is no fill].
+
+    cmap : str
+        Give a CPT. The *r*-value for each sector is used to look-up the
+        sector color. Cannot be used with a rose diagram.
+
+    pen : str
+        Set pen attributes for sector outline or rose plot, e.g.
+        ``pen = '0.5p'``. [Default is no outline]. To change pen used to
+        draw vector (requires ``vectors``) [Default is same as sector
+        outline] use e.g. ``pen = 'v0.5p'``.
+
+    labels : str
+         ``'wlabel,elabel,slabel,nlabel'``. Specify labels for the 0, 90,
+         180, and 270 degree marks. For full-circle plot the default is
+         WEST,EAST,SOUTH,NORTH and for half-circle the default is
+         90W,90E,-,0. A **-** in any entry disables that label
+         (e.g. ``labels = 'W,E,-,N'``). Use ``labels = ''`` to disable
+         all four labels. Note that the :gmt-term:`GMT_LANGUAGE` setting
+         will affect the words used.
+
+    no_scale : bool
+         Do NOT draw the scale length bar (``no_scale = True``).
+         Default plots scale in lower right corner provided ``frame``
+         is used.
+
+    shift : bool
+         Shift sectors so that they are centered on the bin interval
+         (e.g., first sector is centered on 0 degrees).
+
+    vectors : str
+          ``vectors = 'mode_file'``. Plot vectors showing the
+          principal directions given in the *mode_file* file.
+          Alternatively, specify ``vectors`` to compute and plot
+          mean direction. See ``vector_params`` to control the vector
+          attributes. Finally, to instead save the computed mean
+          direction and other statistics, use
+          ``vectors = '+wmode_file'``. The eight items saved to
+          a single record are: *mean_az*, *mean_r*, *mean_resultant*,
+          *max_r*, *scaled_mean_r*, *length_sum*, *n*, *sign@alpha*,
+          where the last term is 0 or 1 depending on whether the mean
+          resultant is significant at the level of confidence set via
+          ``alpha``.
+
+    vector_params : str
+        Used with ``vectors`` to modify vector parameters. For
+        vector heads, append vector head size [Default is 0, i.e., a
+        line]. See :gmt-docs:`rose.html#vector-attributes` for
+        specifying additional attributes. If ``vectors`` is not
+        given and the current plot mode is to draw a windrose diagram
+        then using ``vector_params`` will add vector heads to all
+        individual directions using the supplied attributes.
+
+    alpha : float or str
+        Sets the confidence level used to determine if the mean
+        resultant is significant (i.e., Lord Rayleigh test for
+        uniformity) [``alpha = 0.05``]. Note: The critical
+        values are approximated [Berens, 2009] and requires at
+        least 10 points; the critical resultants are accurate to at
+        least 3 significant digits. For smaller data sets you
+        should consult exact statistical tables.
+
+        Berens, P., 2009, CircStat: A MATLAB Toolbox for Circular
+        Statistics, *J. Stat. Software*, 31(10), 1-21,
+        https://doi.org/10.18637/jss.v031.i10.
+
+    {U}
+    {V}
+    {XY}
+    {c}
+    {p}
+    {t}
+    """
+
+    kwargs = self._preprocess(**kwargs)  # pylint: disable=protected-access
+
+    with Session() as lib:
+        # Choose how data will be passed into the module
+        file_context = lib.virtualfile_from_data(
+            check_kind="vector", data=data, x=length, y=azimuth
+        )
+
+        with file_context as fname:
+            arg_str = " ".join([fname, build_arg_string(kwargs)])
+
+            lib.call_module("rose", arg_str)
diff --git a/pygmt/tests/baseline/test_rose_2d_array_multiple.png.dvc b/pygmt/tests/baseline/test_rose_2d_array_multiple.png.dvc
new file mode 100644
index 00000000000..d2f072536a5
--- /dev/null
+++ b/pygmt/tests/baseline/test_rose_2d_array_multiple.png.dvc
@@ -0,0 +1,4 @@
+outs:
+- md5: 2e440e6db215a619829acb08a0ab55bf
+  size: 41397
+  path: test_rose_2d_array_multiple.png
diff --git a/pygmt/tests/baseline/test_rose_2d_array_single.png.dvc b/pygmt/tests/baseline/test_rose_2d_array_single.png.dvc
new file mode 100644
index 00000000000..886ea8638b3
--- /dev/null
+++ b/pygmt/tests/baseline/test_rose_2d_array_single.png.dvc
@@ -0,0 +1,4 @@
+outs:
+- md5: c64abda4022e5e9fcf5592c5c2f6bb44
+  size: 37742
+  path: test_rose_2d_array_single.png
diff --git a/pygmt/tests/baseline/test_rose_bools.png.dvc b/pygmt/tests/baseline/test_rose_bools.png.dvc
new file mode 100644
index 00000000000..bf45d701a36
--- /dev/null
+++ b/pygmt/tests/baseline/test_rose_bools.png.dvc
@@ -0,0 +1,4 @@
+outs:
+- md5: 2b0ba264c2240a6d2da02466d8848d7a
+  size: 73857
+  path: test_rose_bools.png
diff --git a/pygmt/tests/baseline/test_rose_data_file.png.dvc b/pygmt/tests/baseline/test_rose_data_file.png.dvc
new file mode 100644
index 00000000000..e2788260503
--- /dev/null
+++ b/pygmt/tests/baseline/test_rose_data_file.png.dvc
@@ -0,0 +1,4 @@
+outs:
+- md5: eba3b08913112996eca8b9d01f5709c3
+  size: 40610
+  path: test_rose_data_file.png
diff --git a/pygmt/tests/baseline/test_rose_no_sectors.png.dvc b/pygmt/tests/baseline/test_rose_no_sectors.png.dvc
new file mode 100644
index 00000000000..4b9404913cc
--- /dev/null
+++ b/pygmt/tests/baseline/test_rose_no_sectors.png.dvc
@@ -0,0 +1,4 @@
+outs:
+- md5: 67613d130577c8083e1d505155421205
+  size: 159974
+  path: test_rose_no_sectors.png
diff --git a/pygmt/tests/baseline/test_rose_plot_data_using_cpt.png.dvc b/pygmt/tests/baseline/test_rose_plot_data_using_cpt.png.dvc
new file mode 100644
index 00000000000..36c4d8b7fb7
--- /dev/null
+++ b/pygmt/tests/baseline/test_rose_plot_data_using_cpt.png.dvc
@@ -0,0 +1,4 @@
+outs:
+- md5: 1facbc2eb905197cd6aece39c08c8686
+  size: 44793
+  path: test_rose_plot_data_using_cpt.png
diff --git a/pygmt/tests/baseline/test_rose_plot_with_transparency.png.dvc b/pygmt/tests/baseline/test_rose_plot_with_transparency.png.dvc
new file mode 100644
index 00000000000..aa0e5c5d1fc
--- /dev/null
+++ b/pygmt/tests/baseline/test_rose_plot_with_transparency.png.dvc
@@ -0,0 +1,4 @@
+outs:
+- md5: 11af5a2d9971898c96f48e0a5215d759
+  size: 62490
+  path: test_rose_plot_with_transparency.png
diff --git a/pygmt/tests/test_rose.py b/pygmt/tests/test_rose.py
new file mode 100644
index 00000000000..cef61207f8d
--- /dev/null
+++ b/pygmt/tests/test_rose.py
@@ -0,0 +1,190 @@
+"""
+Tests for rose.
+"""
+import numpy as np
+import pytest
+from pygmt import Figure
+from pygmt.datasets import load_fractures_compilation
+
+
+@pytest.fixture(scope="module", name="data")
+def fixture_data():
+    """
+    Load the sample numpy array data.
+    """
+    return np.array(
+        [[40, 60], [60, 300], [20, 180], [30, 190], [60, 90], [40, 110], [80, 125]]
+    )
+
+
+@pytest.fixture(scope="module", name="data_fractures_compilation")
+def fixture_data_fractures_compilation():
+    """
+    Load the sample fractures compilation dataset.
+    """
+    return load_fractures_compilation()
+
+
+@pytest.mark.mpl_image_compare
+def test_rose_data_file(data_fractures_compilation):
+    """
+    Test supplying data from sample dataset.
+    """
+    fig = Figure()
+    fig.rose(
+        data=data_fractures_compilation,
+        region=[0, 1, 0, 360],
+        sector=15,
+        diameter="5.5c",
+        color="blue",
+        frame=["x0.2g0.2", "y30g30", "+glightgray"],
+        columns=[1, 0],
+        pen="1p",
+        norm="",
+        scale=0.4,
+    )
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+def test_rose_2d_array_single():
+    """
+    Test supplying a 2D numpy array containing a single pair of lengths and
+    directions.
+    """
+    data = np.array([[40, 60]])
+    fig = Figure()
+    fig.rose(
+        data=data,
+        region=[0, 1, 0, 360],
+        sector=10,
+        diameter="5.5c",
+        color="cyan",
+        frame=["x0.2g0.2", "y30g30", "+glightgray"],
+        pen="1p",
+        norm=True,
+        scale=0.4,
+    )
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+def test_rose_2d_array_multiple(data):
+    """
+    Test supplying a 2D numpy array containing a list of lengths and
+    directions.
+    """
+    fig = Figure()
+    fig.rose(
+        data=data,
+        region=[0, 1, 0, 360],
+        sector=10,
+        diameter="5.5c",
+        color="blue",
+        frame=["x0.2g0.2", "y30g30", "+gmoccasin"],
+        pen="1p",
+        norm=True,
+        scale=0.4,
+    )
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+def test_rose_plot_data_using_cpt(data):
+    """
+    Test supplying a 2D numpy array containing a list of lengths and
+    directions.
+
+    Use a cmap to color sectors.
+    """
+    fig = Figure()
+    fig.rose(
+        data=data,
+        region=[0, 1, 0, 360],
+        sector=15,
+        diameter="5.5c",
+        cmap="batlow",
+        frame=["x0.2g0.2", "y30g30", "+gdarkgray"],
+        pen="1p",
+        norm=True,
+        scale=0.4,
+    )
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+def test_rose_plot_with_transparency(data_fractures_compilation):
+    """
+    Test supplying a data file containing a list of fracture lengths and
+    azimuth as digitized from geological maps to the data argument (lengths are
+    stored in the second column, azimuths in the first, specify via columns).
+
+    Use transparency.
+    """
+    fig = Figure()
+    fig.rose(
+        data=data_fractures_compilation,
+        region=[0, 1, 0, 360],
+        sector=15,
+        diameter="5.5c",
+        color="blue",
+        frame=["x0.2g0.2", "y30g30", "+glightgray"],
+        columns=[1, 0],
+        pen="1p",
+        norm=True,
+        scale=0.4,
+        transparency=50,
+    )
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+def test_rose_no_sectors(data_fractures_compilation):
+    """
+    Test supplying a data file containing a list of fracture lengths and
+    azimuth as digitized from geological maps to the data argument (lengths are
+    stored in the second column, azimuths in the first, specify via columns).
+
+    Plot data without defining a sector width, add a title and rename labels.
+    """
+    fig = Figure()
+    fig.rose(
+        data=data_fractures_compilation,
+        region=[0, 500, 0, 360],
+        columns="1,0",
+        diameter="10c",
+        labels="180/0/90/270",
+        frame=["xg100", "yg45", "+t'Windrose diagram'"],
+        pen="1.5p,red3",
+        transparency=40,
+        scale=0.5,
+    )
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+def test_rose_bools(data_fractures_compilation):
+    """
+    Test supplying a data file containing a list of fracture lengths and
+    azimuth as digitized from geological maps to the data argument (lengths are
+    stored in the second column, azimuths in the first, specify via columns).
+
+    Test bools.
+    """
+    fig = Figure()
+    fig.rose(
+        data=data_fractures_compilation,
+        region=[0, 1, 0, 360],
+        sector=10,
+        columns=[1, 0],
+        diameter="10c",
+        frame=["x0.2g0.2", "y30g30", "+glightgray"],
+        color="red3",
+        pen="1p",
+        orientation=False,
+        norm=True,
+        vectors=True,
+        no_scale=True,
+        shift=False,
+    )
+    return fig