diff --git a/docs/getting_started.rst b/docs/getting_started.rst
index 28c411a..abee52b 100644
--- a/docs/getting_started.rst
+++ b/docs/getting_started.rst
@@ -77,6 +77,32 @@ To print out summary information about the result::
 The `~hips.HipsDrawResult` object also gives access to the `~hips.HipsTile`
 objects that were used for drawing the sky image, as well as other things.
 
+Plot using Astropy visualization toolkit
+========================================
+
+Astropy provides a framework for plotting astronomical images with coordinates. It builds on top of Matplotlib and provides functionalities such as image normalization (scaling and stretching), smart histogram plotting, RGB color image creation from separate images. The framework also allows for customization of plotting styles.
+
+The example below is for the FITS format and controls the stretch of the image through normalization. For FITS tiles, the data type is either ``int16`` or ``float32``::
+
+    import matplotlib.pyplot as plt
+    from astropy.visualization.mpl_normalize import simple_norm
+
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+    im = ax.imshow(result.image, origin='lower', norm=norm, cmap='gray')
+    fig.colorbar(im)
+
+.. plot:: plot_fits_astropy.py
+
+RGB tiles can be plotted in much the same way as above, however, it is uncommon to apply an extra stretch in this case. For ``jpg`` and ``png`` tiles, the data type is ``uint8``.
+
+.. note::
+
+    For ``png`` tiles, there are four channel i.e. RGBA. The alpha channel is used for controlling the transparency of the image.
+
+The example provided here is trivial. Astropy provides numerous other features, customizability options, and in-depth examples. Please see their documentation at:
+https://docs.astropy.org/en/stable/visualization
+
 Make a color sky image
 ======================
 
@@ -85,11 +111,10 @@ Making a color sky image works the same as the grayscale image example above,
 except that you get back a 3-dim Numpy array with ``(R, G, B)`` channels for ``jpg``
 or ``(R, G, B, A)`` channels (``A`` is transparency) for ``png``.
 
-Here's an example using ``jpg`` and http://alasky.u-strasbg.fr/Fermi/Color/ :
+Here's an example using ``jpg`` and http://alasky.u-strasbg.fr/Fermi/Color:
 
 .. plot:: plot_jpg.py
 
-
 HiPS data
 =========
 
diff --git a/docs/plot_fits.py b/docs/plot_fits.py
index 392e80d..270e79f 100644
--- a/docs/plot_fits.py
+++ b/docs/plot_fits.py
@@ -11,13 +11,3 @@
 hips_survey = 'CDS/P/DSS2/red'
 result = make_sky_image(geometry=geometry, hips_survey=hips_survey, tile_format='fits')
 result.plot()
-
-# Draw the sky image
-# import matplotlib.pyplot as plt
-# from astropy.visualization.mpl_normalize import simple_norm
-# ax = plt.subplot(projection=geometry.wcs)
-# norm = simple_norm(result.image, 'sqrt', min_percent=1, max_percent=99)
-# ax.imshow(result.image, origin='lower', norm=norm, cmap='gray')
-
-# import matplotlib.pyplot as plt
-# plt.show()
diff --git a/docs/plot_fits_astropy.py b/docs/plot_fits_astropy.py
new file mode 100644
index 0000000..2ffc4d6
--- /dev/null
+++ b/docs/plot_fits_astropy.py
@@ -0,0 +1,25 @@
+"""Basic example how to plot a sky image with the hips package"""
+from astropy.coordinates import SkyCoord
+from hips import WCSGeometry, make_sky_image
+
+# Compute the sky image
+geometry = WCSGeometry.create(
+    skydir=SkyCoord(0, 0, unit='deg', frame='galactic'),
+    width=2000, height=1000, fov="3 deg",
+    coordsys='galactic', projection='AIT',
+)
+hips_survey = 'CDS/P/DSS2/red'
+result = make_sky_image(geometry=geometry, hips_survey=hips_survey, tile_format='fits')
+
+# Draw the sky image
+import matplotlib.pyplot as plt
+from astropy.visualization.mpl_normalize import simple_norm
+
+# Perform normalization.
+norm = simple_norm(result.image, 'sqrt', min_percent=1, max_percent=99)
+
+# Display the image
+fig = plt.figure()
+ax = fig.add_subplot(1, 1, 1)
+im = ax.imshow(result.image, origin='lower', norm=norm, cmap='gray')
+fig.colorbar(im)
diff --git a/hips/draw/ui.py b/hips/draw/ui.py
index 50413d6..b12ee19 100644
--- a/hips/draw/ui.py
+++ b/hips/draw/ui.py
@@ -121,19 +121,27 @@ def write_image(self, filename: str, overwrite: bool = False) -> None:
             image = Image.fromarray(self.image)
             image.save(filename)
 
-    def plot(self) -> None:
+    def plot(self, show_grid: bool = False) -> None:
         """Plot the all sky image and overlay HiPS tile outlines.
 
+        Parameters
+        ----------
+        show_grid : bool
+            Enable grid around HiPS tile boundaries
+
         Uses `astropy.visualization.wcsaxes`.
         """
         import matplotlib.pyplot as plt
-        for tile in self.tiles:
-            corners = tile.meta.skycoord_corners
-            corners = corners.transform_to(self.geometry.celestial_frame)
-            ax = plt.subplot(projection=self.geometry.wcs)
-            opts = dict(color='red', lw=1, )
-            ax.plot(corners.data.lon.deg, corners.data.lat.deg,
-                    transform=ax.get_transform('world'), **opts)
+        ax = plt.subplot(projection=self.geometry.wcs)
+
+        if show_grid:
+            for tile in self.tiles:
+                corners = tile.meta.skycoord_corners
+                corners = corners.transform_to(self.geometry.celestial_frame)
+                opts = dict(color='red', lw=1)
+                ax.plot(corners.data.lon.deg, corners.data.lat.deg,
+                        transform=ax.get_transform('world'), **opts)
+
         ax.imshow(self.image, origin='lower')
 
     def report(self) -> None: