Incorrect label for FluxPoints plots
**Gammapy version**
v0.19 an dev

**Bug description**
When plotting a `FluxPoints` in `dnde` or `e2dnde`, the label is confusing because it gives an interval where it should give the reference energy. 
**Expected behavior**
Give only the reference energy. Possible add it in the y-axis title.

**To Reproduce**
Just look at the flight curve flare tutorial. It is incorrect there.


**Other information**
Any other information you think will be useful for us to fix the issue can go here.

diff --git a/Makefile b/Makefile
index a9a1c456b..4f05a7b17 100644
--- a/Makefile
+++ b/Makefile
@@ -48,7 +48,7 @@ help:
 	@echo ''
 
 clean:
-	rm -rf build dist docs/_build docs/api htmlcov MANIFEST gammapy.egg-info .coverage
+	rm -rf build dist docs/_build docs/api htmlcov MANIFEST gammapy.egg-info .coverage .cache
 	find . -name "*.pyc" -exec rm {} \;
 	find . -name "*.so" -exec rm {} \;
 	find gammapy -name '*.c' -exec rm {} \;
diff --git a/gammapy/catalog/fermi.py b/gammapy/catalog/fermi.py
index 08a83fac9..9e8ec5ec8 100644
--- a/gammapy/catalog/fermi.py
+++ b/gammapy/catalog/fermi.py
@@ -14,7 +14,7 @@ from astropy.modeling.models import Gaussian2D, Disk2D
 from astropy.coordinates import Angle
 from ..utils.scripts import make_path
 from ..utils.energy import EnergyBounds
-from ..utils.table import table_standardise
+from ..utils.table import table_standardise_units_inplace
 from ..image import SkyImage
 from ..image.models import Delta2D, Template2D
 from ..spectrum import FluxPoints, compute_flux_points_dnde
@@ -699,7 +699,7 @@ class SourceCatalog3FGL(SourceCatalog):
 
         with ignore_warnings():  # ignore FITS units warnings
             table = Table.read(filename, hdu='LAT_Point_Source_Catalog')
-        table_standardise(table)
+        table_standardise_units_inplace(table)
 
         source_name_key = 'Source_Name'
         source_name_alias = ('Extended_Source_Name', '0FGL_Name', '1FGL_Name',
@@ -728,7 +728,7 @@ class SourceCatalog1FHL(SourceCatalog):
 
         with ignore_warnings():  # ignore FITS units warnings
             table = Table.read(filename, hdu='LAT_Point_Source_Catalog')
-        table_standardise(table)
+        table_standardise_units_inplace(table)
 
         source_name_key = 'Source_Name'
         source_name_alias = ('ASSOC1', 'ASSOC2', 'ASSOC_TEV', 'ASSOC_GAM')
@@ -755,7 +755,7 @@ class SourceCatalog2FHL(SourceCatalog):
 
         with ignore_warnings():  # ignore FITS units warnings
             table = Table.read(filename, hdu='2FHL Source Catalog')
-        table_standardise(table)
+        table_standardise_units_inplace(table)
 
         source_name_key = 'Source_Name'
         source_name_alias = ('ASSOC', '3FGL_Name', '1FHL_Name', 'TeVCat_Name')
@@ -784,7 +784,7 @@ class SourceCatalog3FHL(SourceCatalog):
 
         with ignore_warnings():  # ignore FITS units warnings
             table = Table.read(filename, hdu='LAT_Point_Source_Catalog')
-        table_standardise(table)
+        table_standardise_units_inplace(table)
 
         source_name_key = 'Source_Name'
         source_name_alias = ('ASSOC1', 'ASSOC2', 'ASSOC_TEV', 'ASSOC_GAM')
diff --git a/gammapy/image/radial_profile.py b/gammapy/image/radial_profile.py
index b7892b0fa..f870fe7cb 100644
--- a/gammapy/image/radial_profile.py
+++ b/gammapy/image/radial_profile.py
@@ -1,6 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import numpy as np
-from astropy.table import QTable
+from astropy.table import Table
 from .core import SkyImage
 
 __all__ = [
@@ -30,7 +30,7 @@ def radial_profile(image, center, radius):
 
     Returns
     -------
-    table : `~astropy.table.QTable`
+    table : `~astropy.table.Table`
 
         Table with the following fields that define the binning:
 
@@ -145,7 +145,7 @@ def _radial_profile_measure(image, labels, index):
     image = image.data
     labels = labels.data
 
-    table = QTable()
+    table = Table()
     table['N_PIX'] = ndimage.sum(np.ones_like(image), labels, index=index)
     table['SUM'] = ndimage.sum(image, labels, index=index)
     table['MEAN'] = table['SUM'] / table['N_PIX']
diff --git a/gammapy/scripts/image_pipe.py b/gammapy/scripts/image_pipe.py
index 16199e2fe..52d9870eb 100755
--- a/gammapy/scripts/image_pipe.py
+++ b/gammapy/scripts/image_pipe.py
@@ -2,7 +2,7 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 import numpy as np
 from astropy.units import Quantity
-from astropy.table import QTable, Table
+from astropy.table import Table
 from astropy.coordinates import Angle
 from astropy.convolution import Tophat2DKernel
 from ..utils.energy import EnergyBounds
@@ -115,7 +115,7 @@ class SingleObsImageMaker(object):
 
         Returns
         -------
-        table : `astropy.table.QTable`
+        table : `astropy.table.Table`
             Two columns: offset in the FOV "theta" and expected counts "npred"
         """
         energy = EnergyBounds.equal_log_spacing(self.energy_band[0].value, self.energy_band[1].value, 100,
@@ -133,7 +133,7 @@ class SingleObsImageMaker(object):
             norm = np.sum(spectrum * energy_band)
             npred /= norm
 
-        table = QTable()
+        table = Table()
         table['theta'] = offset
         table['npred'] = npred
 
diff --git a/gammapy/spectrum/flux_point.py b/gammapy/spectrum/flux_point.py
index b86d5b5bb..85d35aac7 100644
--- a/gammapy/spectrum/flux_point.py
+++ b/gammapy/spectrum/flux_point.py
@@ -9,6 +9,7 @@ from astropy import units as u
 from astropy.io.registry import IORegistryError
 from ..utils.scripts import make_path
 from ..utils.fits import table_from_row_data
+from ..utils.table import table_standardise_units_copy
 
 __all__ = [
     'compute_flux_points_dnde',
@@ -20,29 +21,30 @@ __all__ = [
 
 log = logging.getLogger(__name__)
 
-REQUIRED_COLUMNS = {'dnde': ['e_ref', 'dnde'],
-                    'e2dnde': ['e_ref', 'e2dnde'],
-                    'flux': ['e_min', 'e_max', 'flux'],
-                    'eflux': ['e_min', 'e_max', 'eflux']}
+REQUIRED_COLUMNS = OrderedDict([
+    ('dnde', ['e_ref', 'dnde']),
+    ('e2dnde', ['e_ref', 'e2dnde']),
+    ('flux', ['e_min', 'e_max', 'flux']),
+    ('eflux', ['e_min', 'e_max', 'eflux']),
+])
 
-OPTIONAL_COLUMNS = {'dnde': ['dnde_err', 'dnde_errp', 'dnde_errn',
-                             'dnde_ul', 'is_ul'],
-                    'e2dnde': ['e2dnde_err', 'e2dnde_errp', 'e2dnde_errn',
-                               'e2dnde_ul', 'is_ul'],
-                    'flux': ['flux_err', 'flux_errp', 'flux_errn',
-                             'flux_ul', 'is_ul'],
-                    'eflux': ['eflux_err', 'eflux_errp', 'eflux_errn',
-                              'eflux_ul', 'is_ul']}
+OPTIONAL_COLUMNS = OrderedDict([
+    ('dnde', ['dnde_err', 'dnde_errp', 'dnde_errn', 'dnde_ul', 'is_ul']),
+    ('e2dnde', ['e2dnde_err', 'e2dnde_errp', 'e2dnde_errn', 'e2dnde_ul', 'is_ul']),
+    ('flux', ['flux_err', 'flux_errp', 'flux_errn', 'flux_ul', 'is_ul']),
+    ('eflux', ['eflux_err', 'eflux_errp', 'eflux_errn', 'eflux_ul', 'is_ul']),
+])
 
-DEFAULT_UNIT = {'dnde': u.Unit('cm-2 s-1 TeV-1'),
-                'e2dnde': u.Unit('erg cm-2 s-1'),
-                'flux': u.Unit('cm-2 s-1'),
-                'eflux': u.Unit('erg cm-2 s-1')}
+DEFAULT_UNIT = OrderedDict([
+    ('dnde', u.Unit('cm-2 s-1 TeV-1')),
+    ('e2dnde', u.Unit('erg cm-2 s-1')),
+    ('flux', u.Unit('cm-2 s-1')),
+    ('eflux', u.Unit('erg cm-2 s-1')),
+])
 
 
 class FluxPoints(object):
-    """
-    Flux point object.
+    """Flux points container.
 
     For a complete documentation see :ref:`gadf:flux-points`, for an usage
     example see :ref:`flux-point-computation`.
@@ -58,50 +60,29 @@ class FluxPoints(object):
 
     Examples
     --------
-
     >>> from gammapy.spectrum import FluxPoints
     >>> filename = '$GAMMAPY_EXTRA/test_datasets/spectrum/flux_points/flux_points.fits'
     >>> flux_points = FluxPoints.read(filename)
-    >>> flux_points.show()
-
+    >>> flux_points.plot()
     """
 
     def __init__(self, table):
-        # validate that the table is a valid representation of the given
-        # flux point sed type
-
-        # TODO: this is a temp solution
-        # Make sure we don't have "ph" in units
-        # Should we use a unit equivalency?
-        unit_changes = [
-            ('ph cm-2 s-1', 'cm-2 s-1'),
-            ('ph cm-2 s-1 TeV-1', 'cm-2 s-1 TeV-1'),
-            ('ph cm-2 s-1 MeV-1', 'cm-2 s-1 MeV-1'),
-        ]
-        for colname in table.colnames:
-            for unit_old, unit_new in unit_changes:
-                if (table[colname].unit is not None) and (u.Unit(table[colname].unit) == u.Unit(unit_old)):
-                    table[colname].unit = u.Unit(unit_new)
-
+        table = table_standardise_units_copy(table)
+        # validate that the table is a valid representation
+        # of the given flux point sed type
         self.table = self._validate_table(table)
 
     @property
     def sed_type(self):
-        """
-        Flux points sed type.
+        """SED type (str).
 
-        Returns
-        -------
-        sed_type : str
-            Can be either 'dnde', 'e2dnde', 'flux' or 'eflux'.
+        One of: {'dnde', 'e2dnde', 'flux', 'eflux'}
         """
         return self.table.meta['SED_TYPE']
 
     @staticmethod
     def _guess_sed_type(table):
-        """
-        Guess sed type from table content.
-        """
+        """Guess SED type from table content."""
         valid_sed_types = list(REQUIRED_COLUMNS.keys())
         for sed_type in valid_sed_types:
             required = set(REQUIRED_COLUMNS[sed_type])
@@ -110,17 +91,13 @@ class FluxPoints(object):
 
     @staticmethod
     def _guess_sed_type_from_unit(unit):
-        """
-        Guess sed type from unit.
-        """
+        """Guess SED type from unit."""
         for sed_type, default_unit in DEFAULT_UNIT.items():
             if unit.is_equivalent(default_unit):
                 return sed_type
 
     def _validate_table(self, table):
-        """
-        Validate input flux point table.
-        """
+        """Validate input table."""
         table = Table(table)
         sed_type = table.meta['SED_TYPE']
         required = set(REQUIRED_COLUMNS[sed_type])
@@ -132,101 +109,12 @@ class FluxPoints(object):
         return table
 
     def _get_y_energy_unit(self, y_unit):
-        """
-        Get energy part of the given y unit.
-        """
+        """Get energy part of the given y unit."""
         try:
             return [_ for _ in y_unit.bases if _.physical_type == 'energy'][0]
         except IndexError:
             return u.Unit('TeV')
 
-    def plot(self, ax=None, sed_type=None, energy_unit='TeV', flux_unit=None,
-             energy_power=0, **kwargs):
-        """
-        Plot flux points
-
-        Parameters
-        ----------
-        ax : `~matplotlib.axes.Axes`
-            Axis object to plot on.
-        sed_type : ['dnde', 'flux', 'eflux']
-            Which sed type to plot.
-        energy_unit : str, `~astropy.units.Unit`, optional
-            Unit of the energy axis
-        flux_unit : str, `~astropy.units.Unit`, optional
-            Unit of the flux axis
-        energy_power : int
-            Power of energy to multiply y axis with
-        kwargs : dict
-            Keyword arguments passed to :func:`~matplotlib.pyplot.errorbar`
-
-        Returns
-        -------
-        ax : `~matplotlib.axes.Axes`
-            Axis object
-        """
-        import matplotlib.pyplot as plt
-
-        if ax is None:
-            ax = plt.gca()
-
-        sed_type = sed_type or self.sed_type
-        y_unit = u.Unit(flux_unit or DEFAULT_UNIT[sed_type])
-
-        y = self.table[sed_type].quantity.to(y_unit)
-        x = self.e_ref.to(energy_unit)
-
-        # get errors and ul
-        is_ul = self._is_ul
-        x_err_all = self.get_energy_err(sed_type)
-        y_err_all = self.get_flux_err(sed_type)
-
-        # handle energy power
-        e_unit = self._get_y_energy_unit(y_unit)
-        y_unit = y.unit * e_unit ** energy_power
-        y = (y * np.power(x, energy_power)).to(y_unit)
-
-        y_err, x_err = None, None
-
-        if y_err_all:
-            y_errn = (y_err_all[0] * np.power(x, energy_power)).to(y_unit)
-            y_errp = (y_err_all[1] * np.power(x, energy_power)).to(y_unit)
-            y_err = (y_errn[~is_ul].to(y_unit).value,
-                     y_errp[~is_ul].to(y_unit).value)
-
-        if x_err_all:
-            x_errn, x_errp = x_err_all
-            x_err = (x_errn[~is_ul].to(energy_unit).value,
-                     x_errp[~is_ul].to(energy_unit).value)
-
-        # set flux points plotting defaults
-        kwargs.setdefault('marker', '+')
-        kwargs.setdefault('ls', 'None')
-
-        ebar = ax.errorbar(x[~is_ul].value, y[~is_ul].value, yerr=y_err,
-                           xerr=x_err, **kwargs)
-
-        if is_ul.any():
-            if x_err_all:
-                x_errn, x_errp = x_err_all
-                x_err = (x_errn[is_ul].to(energy_unit).value,
-                         x_errp[is_ul].to(energy_unit).value)
-
-            y_ul = self.table[sed_type + '_ul'].quantity
-            y_ul = (y_ul * np.power(x, energy_power)).to(y_unit)
-
-            y_err = (0.5 * y_ul[is_ul].value, np.zeros_like(y_ul[is_ul].value))
-
-            kwargs.setdefault('c', ebar[0].get_color())
-            ax.errorbar(x[is_ul].value, y_ul[is_ul].value, xerr=x_err, yerr=y_err,
-                        uplims=True, **kwargs)
-
-        ax.set_xscale('log', nonposx='clip')
-        ax.set_yscale('log', nonposy='clip')
-        ax.set_xlabel('Energy ({})'.format(energy_unit))
-        ax.set_ylabel('{} ({})'.format(self.sed_type, y_unit))
-        return ax
-
     def get_energy_err(self, sed_type=None):
         """Compute energy error for given sed type"""
         # TODO: sed_type is not used
@@ -246,7 +134,7 @@ class FluxPoints(object):
         if sed_type is None:
             sed_type = self.sed_type
         try:
-            # assymmetric error
+            # asymmetric error
             y_errn = self.table[sed_type + '_errn'].quantity
             y_errp = self.table[sed_type + '_errp'].quantity
             y_err = (y_errn, y_errp)
@@ -267,29 +155,6 @@ class FluxPoints(object):
         except KeyError:
             return np.isnan(self.table[self.sed_type])
 
-    def peek(self, figsize=(8, 5), **kwargs):
-        """
-        Show flux points.
-
-        Parameters
-        ----------
-        figsize : tuple
-            Figure size
-        kwargs : dict
-            Keyword arguments passed to `FluxPoints.plot()`.
-
-        Returns
-        -------
-        ax : `~matplotlib.axes.Axes`
-            Plotting axes object.
-        """
-        import matplotlib.pyplot as plt
-
-        fig = plt.figure(figsize=figsize)
-        ax = fig.add_subplot(111)
-        self.plot(ax=ax, **kwargs)
-        return ax
-
     def __str__(self):
         """
         String representation of the flux points class.
@@ -306,8 +171,7 @@ class FluxPoints(object):
 
     @classmethod
     def read(cls, filename, **kwargs):
-        """
-        Read flux points.
+        """Read flux points.
 
         Parameters
         ----------
@@ -331,8 +195,7 @@ class FluxPoints(object):
         return cls(table=table)
 
     def write(self, filename, **kwargs):
-        """
-        Write flux points.
+        """Write flux points.
 
         Parameters
         ----------
@@ -348,11 +211,9 @@ class FluxPoints(object):
             kwargs.setdefault('format', 'ascii.ecsv')
             self.table.write(str(filename), **kwargs)
 
-    # TODO: handle with Energy or EnergyBounds classes?
     @property
     def e_ref(self):
-        """
-        Reference energy.
+        """Reference energy.
 
         Defined by `e_ref` column in `FluxPoints.table` or computed as log
         center, if `e_min` and `e_max` columns are present in `FluxPoints.table`.
@@ -365,14 +226,11 @@ class FluxPoints(object):
         try:
             return self.table['e_ref'].quantity
         except KeyError:
-            e_ref = np.sqrt(self.e_min * self.e_max)
-            return e_ref
+            return np.sqrt(self.e_min * self.e_max)
 
-    # TODO: handle with Energy or EnergyBounds classes?
     @property
     def e_min(self):
-        """
-        Lower bound of energy bin.
+        """Lower bound of energy bin.
 
         Defined by `e_min` column in `FluxPoints.table`.
 
@@ -383,11 +241,9 @@ class FluxPoints(object):
         """
         return self.table['e_min'].quantity
 
-    # TODO: handle with Energy or EnergyBounds classes?
     @property
     def e_max(self):
-        """
-        Upper bound of energy bin.
+        """Upper bound of energy bin.
 
         Defined by `e_max` column in `FluxPoints.table`.
 
@@ -399,24 +255,23 @@ class FluxPoints(object):
         return self.table['e_max'].quantity
 
     def drop_ul(self):
-        """
-        Drop upper limit flux points.
-
-        Examples
-        --------
-
-        from gammapy.spectrum import FluxPoints
-
-        filename = '$GAMMAPY_EXTRA/test_datasets/spectrum/flux_points/flux_points.fits'
-        flux_points = FluxPoints.read(filename)
-
-        print(flux_points)
-        print(flux_points.drop_ul())
+        """Drop upper limit flux points.
 
         Returns
         -------
         flux_points : `FluxPoints`
             Flux points with upper limit points removed.
+
+        Examples
+        --------
+
+        ::
+
+            from gammapy.spectrum import FluxPoints
+            filename = '$GAMMAPY_EXTRA/test_datasets/spectrum/flux_points/flux_points.fits'
+            flux_points = FluxPoints.read(filename)
+            print(flux_points)
+            print(flux_points.drop_ul())
         """
         table = self.table.copy()
         table_drop_ul = table[~self._is_ul]
@@ -424,16 +279,14 @@ class FluxPoints(object):
 
     @classmethod
     def stack(cls, flux_points):
-        """
-        Create a new `FluxPoints` object by stacking a list of existing
-        flux points.
+        """Create flux points by stacking list of flux points.
 
         The first `FluxPoints` object in the list is taken as a reference to infer
         column names and units for the stacked object.
 
         Parameters
         ----------
-        flux_points : list of `FluxPoints` objects
+        flux_points : list of `FluxPoints`
             List of flux points to stack.
 
         Returns
@@ -453,8 +306,95 @@ class FluxPoints(object):
             tables.append(table[reference.colnames])
         table_stacked = vstack(tables)
         table_stacked.meta['SED_TYPE'] = reference.meta['SED_TYPE']
+
         return cls(table_stacked)
 
+    def plot(self, ax=None, sed_type=None, energy_unit='TeV', flux_unit=None,
+             energy_power=0, **kwargs):
+        """Plot flux points.
+
+        Parameters
+        ----------
+        ax : `~matplotlib.axes.Axes`
+            Axis object to plot on.
+        sed_type : ['dnde', 'flux', 'eflux']
+            Which sed type to plot.
+        energy_unit : str, `~astropy.units.Unit`, optional
+            Unit of the energy axis
+        flux_unit : str, `~astropy.units.Unit`, optional
+            Unit of the flux axis
+        energy_power : int
+            Power of energy to multiply y axis with
+        kwargs : dict
+            Keyword arguments passed to :func:`~matplotlib.pyplot.errorbar`
+
+        Returns
+        -------
+        ax : `~matplotlib.axes.Axes`
+            Axis object
+        """
+        import matplotlib.pyplot as plt
+
+        if ax is None:
+            ax = plt.gca()
+
+        sed_type = sed_type or self.sed_type
+        y_unit = u.Unit(flux_unit or DEFAULT_UNIT[sed_type])
+
+        y = self.table[sed_type].quantity.to(y_unit)
+        x = self.e_ref.to(energy_unit)
+
+        # get errors and ul
+        is_ul = self._is_ul
+        x_err_all = self.get_energy_err(sed_type)
+        y_err_all = self.get_flux_err(sed_type)
+
+        # handle energy power
+        e_unit = self._get_y_energy_unit(y_unit)
+        y_unit = y.unit * e_unit ** energy_power
+        y = (y * np.power(x, energy_power)).to(y_unit)
+
+        y_err, x_err = None, None
+
+        if y_err_all:
+            y_errn = (y_err_all[0] * np.power(x, energy_power)).to(y_unit)
+            y_errp = (y_err_all[1] * np.power(x, energy_power)).to(y_unit)
+            y_err = (y_errn[~is_ul].to(y_unit).value,
+                     y_errp[~is_ul].to(y_unit).value)
+
+        if x_err_all:
+            x_errn, x_errp = x_err_all
+            x_err = (x_errn[~is_ul].to(energy_unit).value,
+                     x_errp[~is_ul].to(energy_unit).value)
+
+        # set flux points plotting defaults
+        kwargs.setdefault('marker', '+')
+        kwargs.setdefault('ls', 'None')
+
+        ebar = ax.errorbar(x[~is_ul].value, y[~is_ul].value, yerr=y_err,
+                           xerr=x_err, **kwargs)
+
+        if is_ul.any():
+            if x_err_all:
+                x_errn, x_errp = x_err_all
+                x_err = (x_errn[is_ul].to(energy_unit).value,
+                         x_errp[is_ul].to(energy_unit).value)
+
+            y_ul = self.table[sed_type + '_ul'].quantity
+            y_ul = (y_ul * np.power(x, energy_power)).to(y_unit)
+
+            y_err = (0.5 * y_ul[is_ul].value, np.zeros_like(y_ul[is_ul].value))
+
+            kwargs.setdefault('c', ebar[0].get_color())
+            ax.errorbar(x[is_ul].value, y_ul[is_ul].value, xerr=x_err, yerr=y_err,
+                        uplims=True, **kwargs)
+
+        ax.set_xscale('log', nonposx='clip')
+        ax.set_yscale('log', nonposy='clip')
+        ax.set_xlabel('Energy ({})'.format(energy_unit))
+        ax.set_ylabel('{} ({})'.format(self.sed_type, y_unit))
+        return ax
+
 
 def compute_flux_points_dnde(flux_points, model, method='lafferty'):
     """
diff --git a/gammapy/spectrum/tests/test_flux_point.py b/gammapy/spectrum/tests/test_flux_point.py
index a2de83ba2..2b428d894 100644
--- a/gammapy/spectrum/tests/test_flux_point.py
+++ b/gammapy/spectrum/tests/test_flux_point.py
@@ -246,12 +246,11 @@ class TestSEDLikelihoodProfile:
         self.sed = SEDLikelihoodProfile.read(filename)
 
     def test_basics(self):
-        # print(self.sed)
         assert 'SEDLikelihoodProfile' in str(self.sed)
 
     @requires_dependency('matplotlib')
     def test_plot(self):
-        ax = self.sed.plot()
+        self.sed.plot()
 
 
 @pytest.fixture(params=FLUX_POINTS_FILES, scope='session')
@@ -263,12 +262,10 @@ def flux_points(request):
 @requires_dependency('yaml')
 @requires_data('gammapy-extra')
 class TestFluxPoints:
+
     @requires_dependency('matplotlib')
     def test_plot(self, flux_points):
-        import matplotlib.pyplot as plt
-        ax = plt.gca()
-        flux_points.plot(ax=ax)
-        flux_points.peek()
+        flux_points.plot()
 
     def test_info(self, flux_points):
         info = str(flux_points)
diff --git a/gammapy/utils/fits.py b/gammapy/utils/fits.py
index 3c50865ef..b4eda36b3 100644
--- a/gammapy/utils/fits.py
+++ b/gammapy/utils/fits.py
@@ -222,6 +222,9 @@ def fits_header_to_meta_dict(header):
     return meta
 
 
+# TODO: remove type = 'qtable' to avoid issues?
+# see https://github.com/astropy/astropy/issues/6098
+# see https://github.com/gammapy/gammapy/issues/980
 def table_from_row_data(rows, type='qtable', **kwargs):
     """Helper function to create table objects from row data.
 
diff --git a/gammapy/utils/table.py b/gammapy/utils/table.py
index 38464e77a..e55767038 100644
--- a/gammapy/utils/table.py
+++ b/gammapy/utils/table.py
@@ -2,16 +2,42 @@
 """Table helper utilities.
 """
 from __future__ import absolute_import, division, print_function, unicode_literals
+from astropy.table import Table
 from .units import standardise_unit
 
 __all__ = [
-    'table_standardise',
+    'table_standardise_units_copy',
+    'table_standardise_units_inplace',
 ]
 
 
-def table_standardise(table):
-    """TODO
+def table_standardise_units_copy(table):
+    """Standardise units for all columns in a table in a copy.
+    
+    Calls `~gammapy.utils.units.standardise_unit`.
+
+    Parameters
+    ----------
+    table : `~astropy.table.Table` or `~astropy.table.QTable`
+        Input table (won't be modified)
+
+    Returns
+    -------
+    table : `~astropy.table.Table`
+        Copy of the input table with standardised column units
+    """
+    # Note: we could add an `inplace` option (or variant of this function)
+    # for `Table`, but not for `QTable`.
+    # See https://github.com/astropy/astropy/issues/6098
+    table = Table(table)
+    return table_standardise_units_inplace(table)
+
+
+def table_standardise_units_inplace(table):
+    """Standardise units for all columns in a table in place.
     """
     for column in table.columns.values():
         if column.unit:
             column.unit = standardise_unit(column.unit)
+
+    return table
diff --git a/gammapy/utils/tests/test_table.py b/gammapy/utils/tests/test_table.py
index aa2677d09..01bac77b0 100644
--- a/gammapy/utils/tests/test_table.py
+++ b/gammapy/utils/tests/test_table.py
@@ -1,18 +1,20 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import absolute_import, division, print_function, unicode_literals
+from astropy.tests.helper import pytest
 import astropy.units as u
-from astropy.table import Table, Column
-from ..table import table_standardise
+from astropy.table import Table, QTable, Column
+from ..table import table_standardise_units_copy
 
 
-def test_table_standardise():
-    table = Table()
+@pytest.mark.parametrize('table_class', [Table, QTable])
+def test_table_standardise_units(table_class):
+    table = table_class()
     table['a'] = Column([1], unit='ph cm-2 s-1')
     table['b'] = Column([1], unit='ct cm-2 s-1')
     table['c'] = Column([1], unit='cm-2 s-1')
     table['d'] = Column([1])
 
-    table_standardise(table)
+    table = table_standardise_units_copy(table)
 
     assert table['a'].unit == u.Unit('cm-2 s-1')
     assert table['b'].unit == u.Unit('cm-2 s-1')
['Makefile', 'gammapy/scripts/image_pipe.py', 'gammapy/utils/fits.py', 'gammapy/spectrum/flux_point.py', 'gammapy/catalog/fermi.py', 'gammapy/utils/tests/test_table.py', 'gammapy/image/radial_profile.py', 'gammapy/utils/table.py', 'gammapy/spectrum/tests/test_flux_point.py']