diff --git a/pygmt/src/binstats.py b/pygmt/src/binstats.py
index 68f8996b39f..4061c56785b 100644
--- a/pygmt/src/binstats.py
+++ b/pygmt/src/binstats.py
@@ -21,6 +21,7 @@
     N="normalize",
     R="region",
     S="search_radius",
+    T="tiling",
     V="verbose",
     W="weight",
     a="aspatial",
@@ -87,6 +88,21 @@ def binstats(data, **kwargs):
         Set the *search_radius* that determines which data points are
         considered close to a node. Append the distance unit.
         Not compatible with ``tiling``.
+    tiling : str
+        **h**\|\ **r**.
+        Instead of circular, possibly overlapping areas, select
+        non-overlapping tiling. Choose between **r**\ ectangular and
+        **h**\ exagonal binning.
+        For **r**, set bin sizes via ``spacing`` and we write
+        the computed statistics to the grid file named in ``outgrid``.
+        For **h**, we write a table with the centers of the hexagons and
+        the computed statistics to standard output (or to the file named
+        in ``outgrid``). Here, the ``spacing`` setting is expected to
+        set the y-increment only and we compute the x-increment given
+        the geometry. Because the spacing between hexagon centers in
+        x and y directions have a ratio of :math:`\sqrt{{ 3 }}`, we will
+        automatically adjust *xmax* in ``region`` to fit a whole number
+        of hexagons. **Note**: Hexagonal tiling requires Cartesian data.
     weight : str
         Input data have an extra column containing observation point weight.
         If weights are given then weighted statistical quantities will be