diff --git a/lib/matplotlib/colors.py b/lib/matplotlib/colors.py
index c5db6117f1bc..a9a0e2718092 100644
--- a/lib/matplotlib/colors.py
+++ b/lib/matplotlib/colors.py
@@ -1594,55 +1594,14 @@ def linthresh(self, value):
         self._scale.linthresh = value
 
 
+@make_norm_from_scale(
+    scale.PowerScale,
+    init=lambda gamma, vmin=None, vmax=None, clip=False: None)
 class PowerNorm(Normalize):
     """
     Linearly map a given value to the 0-1 range and then apply
     a power-law normalization over that range.
     """
-    def __init__(self, gamma, vmin=None, vmax=None, clip=False):
-        super().__init__(vmin, vmax, clip)
-        self.gamma = gamma
-
-    def __call__(self, value, clip=None):
-        if clip is None:
-            clip = self.clip
-
-        result, is_scalar = self.process_value(value)
-
-        self.autoscale_None(result)
-        gamma = self.gamma
-        vmin, vmax = self.vmin, self.vmax
-        if vmin > vmax:
-            raise ValueError("minvalue must be less than or equal to maxvalue")
-        elif vmin == vmax:
-            result.fill(0)
-        else:
-            if clip:
-                mask = np.ma.getmask(result)
-                result = np.ma.array(np.clip(result.filled(vmax), vmin, vmax),
-                                     mask=mask)
-            resdat = result.data
-            resdat -= vmin
-            resdat[resdat < 0] = 0
-            np.power(resdat, gamma, resdat)
-            resdat /= (vmax - vmin) ** gamma
-
-            result = np.ma.array(resdat, mask=result.mask, copy=False)
-        if is_scalar:
-            result = result[0]
-        return result
-
-    def inverse(self, value):
-        if not self.scaled():
-            raise ValueError("Not invertible until scaled")
-        gamma = self.gamma
-        vmin, vmax = self.vmin, self.vmax
-
-        if np.iterable(value):
-            val = np.ma.asarray(value)
-            return np.ma.power(val, 1. / gamma) * (vmax - vmin) + vmin
-        else:
-            return pow(value, 1. / gamma) * (vmax - vmin) + vmin
 
 
 class BoundaryNorm(Normalize):
diff --git a/lib/matplotlib/scale.py b/lib/matplotlib/scale.py
index 62dc5192a30d..2fea028cd4f1 100644
--- a/lib/matplotlib/scale.py
+++ b/lib/matplotlib/scale.py
@@ -458,6 +458,79 @@ def get_transform(self):
         return self._transform
 
 
+class PowerTransform(Transform):
+    input_dims = output_dims = 1
+
+    def __init__(self, gamma, nonpositive='clip'):
+        super().__init__()
+
+        self.gamma = gamma
+        self._clip = _api.check_getitem(
+            {"clip": True, "mask": False}, nonpositive=nonpositive)
+
+    def __str__(self):
+        return "{}(gamma={}, nonpositive={!r})".format(
+            type(self).__name__, self.gamma, "clip" if self._clip else "mask")
+
+    def transform_non_affine(self, a):
+        with np.errstate(divide="ignore", invalid="ignore"):
+            out = np.power(a, self.gamma)
+            if self._clip:
+                out[a <= 0] = 0
+            return out
+
+    def inverted(self):
+        return InvertedPowerTransform(self.gamma)
+
+
+class InvertedPowerTransform(Transform):
+    input_dims = output_dims = 1
+
+    def __init__(self, gamma):
+        super().__init__()
+        self.gamma = gamma
+
+    def transform_non_affine(self, a):
+        if self.gamma == 0:
+            return np.inf
+        else:
+            return np.power(a, 1./self.gamma)
+
+
+class PowerScale(ScaleBase):
+
+    name = 'power'
+
+    def __init__(self, axis, gamma=0.5):
+        self._transform = PowerTransform(gamma)
+
+    gamma = property(lambda self: self._transform.gamma)
+
+    def get_transform(self):
+
+        return self._transform
+
+    def set_default_locators_and_formatters(self, axis):
+        # docstring inherited
+        axis.set_major_locator(AutoLocator())
+        axis.set_major_formatter(ScalarFormatter())
+        axis.set_minor_formatter(NullFormatter())
+        # update the minor locator for x and y axis based on rcParams
+        if (axis.axis_name == 'x' and mpl.rcParams['xtick.minor.visible'] or
+                axis.axis_name == 'y' and mpl.rcParams['ytick.minor.visible']):
+            axis.set_minor_locator(AutoMinorLocator())
+        else:
+            axis.set_minor_locator(NullLocator())
+
+    def limit_range_for_scale(self, vmin, vmax, minpos):
+        """Limit the domain to positive values."""
+        if not np.isfinite(minpos):
+            minpos = 1e-300
+
+        return (minpos if vmin <= 0 else vmin,
+                minpos if vmax <= 0 else vmax)
+
+
 class LogitTransform(Transform):
     input_dims = output_dims = 1