How to change ColorFunction after plotting

Inspection

Let's take a look at what is being returned from the plot:

plot = DensityPlot[Exp[-((x-5)^2 + (y-5)^2)/5^2], {x,0,10}, {y,0,10}, PlotRange -> {0,1}]

The output mainly consists of a big GraphicsComplex:

In[2]:= Head[Plot[[1]]]
Out[2]= GraphicsComplex

with color specified by VertexColors:

In[3]:= Short[Options[plot[[1]]]]
Out[3]= {VertexColors->{<<1>>}}

So in theory, by replacing VertexColors with other color values, you can change the colors without reevaluation.

Good News

By default, DensityPlot uses color functions on z values and it is scaled (ColorFunctionScaling->True), which makes our task a bit easier since scaled z-values between 0 and 1 will be applied to color functions.

Bad News

We need z-values because we need to apply color functions on them. But, the color function is already evaluated, and only RGB values are stored in VertexColors. Thankfully, the default color function is 1-1 function, but there is no such guarantee in general. Also, the inverse function of the color functions are not trivial (not that hard, just will take some time to compute--beating the purpose of not re-evaluating).

Solution

A bit of cheating. But let's set up the graphics a bit so that it actually contains useful information. Instead of the default color function, let's use GrayLevel.

plot = DensityPlot[Exp[-((x-5)^2 + (y-5)^2)/5^2], {x,0,10}, {y,0,10},
  PlotRange -> {0,1}, ColorFunction -> GrayLevel]

If you take a look at the vertex colors again, now it is in fact the same value as the scaled z-values of the function (and nicely packed too).

In[5]:= Short[Options[plot[[1]]], 4]
Out[5]= {VertexColors->{0.,0.0475548,0.0989165,0.150418,0.197553,0.23559,<<238>>,0.0113745,0.0113745,0.0113745,0.0113745,0.0113745}}

From here, now it is just a matter of applying a color function to the list of vertex colors--which is indeed scaled z-values--and replace it with the new values.

The position of VertexColors within output can be inspected by the following:

In[6]:= Position[Plot, HoldPattern[VertexColors -> _]]
Out[6] = {{1, 3}}

Usually, the positions in plot outputs are pretty stable unless you change options that modifies geometry (such as Exclusions), and there is a slight worry that Position may cause unpacking of nicely packed plot outputs (read this), so let's stick to this value for a while (of course, actual replace should be {1, 3, 2}, the second argument of the rule). A simple test:

ReplacePart[plot, {1,3,2} -> (ColorData["Rainbow"] /@ plot[[1, 3, 2]])]

It works nicely. Now with Manipulate:

Manipulate[
 ReplacePart[plot, {1,3,2} -> (ColorData[gradient] /@ plot[[1, 3, 2]])],
 {gradient, ColorData["Gradients"]}]

Problems

If the goal is to completely avoid Kernel evaluation, there is no easy solution--if not impossible, since ColorData will rely on the Kernel evaluation anyway.

Calling ColorData (and any other color functions, such as Blend or Hue for that matter) manually like this is not efficient in two ways;

First, internal plot functions do a lot of optimization for calling color data. Essentially it is getting a treatment like CompliedFunction, but more specialized. By manually calling them, you are not getting much speed up.

Secondly, the output form is a list of RGBColor[...] which makes the whole result unpacked. Not efficient memory-wise.

There is a way to achieve maximum efficiency by turning ColorData into a linear interpolation function of triples and using CompiledFunction or InterpolatingFunction, but I frankly don't think that it is worth...

One solution that works quite well and is fast (interpolation would've been slow), is to use Nearest on a fine grid. The advantage is that you don't have to re-evaluate your plot, or modify it prior to generating it like Yu-Sung's answer does. Here's an example:

(* Create a rule list of original RGB colors to value and a nearest function *)
list = List @@ ColorData["LakeColors"][#] -> # & /@ Range[0, 1, 0.001];
nf = Nearest[list];

(* function to convert original RGB triplets to those of chosen color data *)
convert[s_String][x_List] := List @@ ColorData[s][First@nf[x]]

(* Apply to plot *)
plot /. Rule[VertexColors, x_List] :> Rule[VertexColors, (convert["Temperature"] /@ x)]

The basic idea is to find the inverse function of the ColorFunction that your original plot uses, to reconstruct the scaled function values. Then apply the new ColorFunction to those values.

Here is how this would look for the special case when the original plot uses Hue:

colorfunction1 = Hue;
colorfunction2 = ColorData["Temperature"]; 
plot = 
 DensityPlot[
  Exp[-((x - 5)^2 + (y - 5)^2)/5^2], {x, 0, 10}, {y, 0, 10}, 
  PlotRange -> {0, 1}, ColorFunction -> colorfunction1]

Now I replace the Hue function:

p = plot /. 
  r_colorfunction1 :> 
   colorfunction2[Quiet[InverseFunction[colorfunction1][r]]] 

So in principle, this works. However, Hue is special because it is invertible easily. For the default color scheme ("LakeColors") that's not so easy. InverseFunction seems to be unavailable for the blends appearing there.

So one would have to do the inversion by hand, probably using FindMinimum, applied to the difference between RGBColors in the VertexColors list and the RGBColors returned by the original color function. It's a FindMinimum problem because numerically the inversion may not work.

Edit

What I did with the Hue example was cheating in the same way that Yu-Sung Chang's solution based on GrayLevel was cheating, because the default color functions are more complicated. So here is something that works for the blends in the ColorData list, including the default LakeColors:

invertColor[rgb_, blend_, x_?NumericQ] := (#.#) &[
   Apply[List, rgb] - Apply[List, blend[x]]];
changeColorFunction[pl_, colfn1_, colfn2_] := 
 pl /. HoldPattern[
    VertexColors -> l_] :> (VertexColors -> Map[colfn2[
        Clip[
         x /. Last@
           Quiet@FindMinimum[invertColor[#, colfn1, x], {x, 0}], {0, 
          1}]] &, l])

The arguments in changeColorFunction[pl_, colfn1_, colfn2_] are the plot that you want to change, the original ColorFunction and the new one (both are assumed to have RGBColor values).

Here I'm applying it to a plot with default colors:

colorfunction2 = ColorData["Temperature"];
colorfunction1 = ColorData["LakeColors"];
plot = DensityPlot[
  Exp[-((x - 5)^2 + (y - 5)^2)/5^2], {x, 0, 10}, {y, 0, 10}, 
  PlotRange -> {0, 1}];
changeColorFunction[plot, colorfunction1, colorfunction2]

Another example with the same default plot:

changeColorFunction[plot, colorfunction1, ColorData["FallColors"]]

Probably it's worth reiterating that the inversion of the color function is a little tricky - but it seems to work well for the default colors.