Ways to plot interpolating functions more economically?

Question

I am interpolating a lot of data over geographic coordinates obtained via GPS, and frequently I need to plot these interpolated functions together. For example, I may have an interpolated elevation map that looks like this:

f0 = InterpolatingFunction[{{x0min, x0max},{y0min, y0max}}, <>]

and several other interpolated functions of other data sets that stretch over subregions of the rectangle {{x0min, x0max}, {y0min, y0max}}. So, say,

f1 = InterpolatingFunction[{{x1min, x1max},{y1min, y1max}}, <>]

where {{x1min, x1max}, {y1min, y1max}} lies entirely within {{x0min, x0max}, {y0min, y0max}} and so on.

I need to plot these interpolated functions A LOT and in attempts to change the visualisations of my data sets, I need an efficient way to extract their limits and pass them on to plotting functions. For this I have written the following function:

intLimits[intF_ /; (Head@intF === InterpolatingFunction), x_: x, y_: y] := 
     Module[{ArgList},
        ArgList = First@intF;
        Which[Length@ArgList < 2,
              ArgList~Prepend~x,
              Length@ArgList == 2,
              {First@ArgList~Prepend~x, Last@ArgList~Prepend~y}
        ]
      ]

which I use with plotting functions like so

f0plot= Plot3D[f0[x,y], Evaluate@First@intLimits[f0], Evaluate@Last@intLimits[f0]];
f1plot= Plot3D[f1[x,y], Evaluate@First@intLimits[f1], Evaluate@Last@intLimits[f1]];

This works well if used to plot a list of interpolation functions (replacing f0 with # and mapping at the list) but it seems a bit ugly and the length of adding Evaluate@First@intLimits doesn't help with the readability of my notebook although seems necessary given the HoldAll Attribute of all plotting functions.

So my question is: is there a more clever way to define the function intLimits or a better way in general to be able to pass on the range of variables of interpolated functions to plotting functions?

Answer 1

Ok, here is one possibility out of many.The following is a function generator:

ClearAll[generate3DPlottingFunction]
generate3DPlottingFunction[fname_Symbol, builtin_Symbol] :=
  Block[{x, y},
     ClearAll[fname];
     fname::no2D = 
         "A function `1` does not represent a 2-dimensional interpolation";
     fname[f_InterpolatingFunction, x_Symbol: x, y_Symbol: y, opts : OptionsPattern[]] :=
       With[{limits = intLimits[f, x, y]},
          With[{fst = First@limits, sec = Last@limits},
             builtin[f[x, y], fst, sec, Evaluate@ FilterRules[{opts}, Options[builtin]]]
          ] /; Length[limits] == 2
       ];
     fname[f_InterpolatingFunction, x_Symbol: x, y_Symbol: y] :=
        (Message[fname::no2D, f]; $Failed);
 fname[___] := $Failed
 ]

Now you can generate the custom functions:

generate3DPlottingFunction[plot3DInterpolated, Plot3D]
generate3DPlottingFunction[contourPlotInterpolated, ContourPlot]

Here is some test function:

testF = 
 Interpolation[
   Flatten[
      Outer[{#1, #2, Sin[#1*#2]} &, Range[0, 2 Pi, 0.1],Range[0, 2 Pi, 0.1]], 
      1
   ]]

Now,

plot3DInterpolated[testF]

contourPlotInterpolated[testF]

Answer 2

I find that the implementation of intLimits[] in the OP is a tad too cumbersome for my taste. It is apparently not too well known that InterpolatingFunction[] objects in Mathematica are set up such that one can extract a number of useful properties from them (I see that rcollyer has alluded to them in a comment). To wit, here's a vastly simpler way to write intLimits[]:

intLimits[intF_InterpolatingFunction, x_Symbol: x, y_Symbol: y] :=
   Join[{{x}, {y}}, intF["Domain"], 2]

One can indeed pass on implementing intLimits[] as a separate function, and just implement a plotting function directly:

interpolantPlot[intFun_InterpolatingFunction, plotFun_Symbol, opts___] := 
     plotFun[intFun[\[FormalX], \[FormalY]], ##, 
             Evaluate[Sequence @@ FilterRules[{opts}, Options[plotFun]]]] & @@
     Join[{{\[FormalX]}, {\[FormalY]}}, intFun["Domain"], 2] /; 
     MatchQ[plotFun, Plot3D | ContourPlot | DensityPlot]

Test it out:

testF = ListInterpolation[Table[Sin[x + Sin[y]], {x, 0, 2 Pi, Pi/30}, {y, 0, 2 Pi, Pi/30}],
                          {{0, 2 Pi}, {0, 2 Pi}}];

interpolantPlot[testF, Plot3D, BoundaryStyle -> Gray, Mesh -> False]

interpolantPlot[testF, ContourPlot, ColorFunction -> "ThermometerColors",
                ContourStyle -> Directive[Gray, Dashed]]