diff --git a/src/Plots.jl b/src/Plots.jl
index 7a97819a..acda8aab 100644
--- a/src/Plots.jl
+++ b/src/Plots.jl
@@ -175,6 +175,7 @@ end
 @shorthands vline
 @shorthands ohlc
 @shorthands contour
+@shorthands contourf
 @shorthands contour3d
 @shorthands surface
 @shorthands wireframe
diff --git a/src/args.jl b/src/args.jl
index 82dd16b5..feb0f32c 100644
--- a/src/args.jl
+++ b/src/args.jl
@@ -79,7 +79,7 @@ add_non_underscore_aliases!(_typeAliases)
 
 like_histogram(seriestype::Symbol) = seriestype in (:histogram, :density)
 like_line(seriestype::Symbol)      = seriestype in (:line, :path, :steppre, :steppost)
-like_surface(seriestype::Symbol)   = seriestype in (:contour, :contour3d, :heatmap, :surface, :wireframe, :image)
+like_surface(seriestype::Symbol)   = seriestype in (:contour, :contourf, :contour3d, :heatmap, :surface, :wireframe, :image)
 
 is3d(seriestype::Symbol) = seriestype in _3dTypes
 is3d(series::Series) = is3d(series.d)
diff --git a/src/backends/pyplot.jl b/src/backends/pyplot.jl
index 339db4a6..ec35af47 100644
--- a/src/backends/pyplot.jl
+++ b/src/backends/pyplot.jl
@@ -423,7 +423,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
             if d[:line_z] == nothing
                 handle = ax[:plot](xyargs...;
                     label = d[:label],
-                    zorder = plt.n,
+                    zorder = d[:series_plotindex],
                     color = py_linecolor(d),
                     linewidth = py_dpi_scale(plt, d[:linewidth]),
                     linestyle = py_linestyle(st, d[:linestyle]),
@@ -504,7 +504,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
     #     end
     #     handle = ax[isvertical(d) ? :bar : :barh](x, y;
     #         label = d[:label],
-    #         zorder = plt.n,
+    #         zorder = d[:series_plotindex],
     #         color = py_fillcolor(d),
     #         edgecolor = py_linecolor(d),
     #         linewidth = d[:linewidth],
@@ -518,7 +518,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
     #     extrakw[isvertical(d) ? :width : :height] = 0.0
     #     handle = ax[isvertical(d) ? :bar : :barh](x, y;
     #         label = d[:label],
-    #         zorder = plt.n,
+    #         zorder = d[:series_plotindex],
     #         color = py_linecolor(d),
     #         edgecolor = py_linecolor(d),
     #         linewidth = d[:linewidth],
@@ -552,7 +552,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
         end
         handle = ax[:scatter](xyargs...;
             label = d[:label],
-            zorder = plt.n + 0.5,
+            zorder = d[:series_plotindex] + 0.5,
             marker = py_marker(d[:markershape]),
             s = py_dpi_scale(plt, d[:markersize] .^ 2),
             edgecolors = py_markerstrokecolor(d),
@@ -565,7 +565,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
     # if st == :histogram
     #     handle = ax[:hist](y;
     #         label = d[:label],
-    #         zorder = plt.n,
+    #         zorder = d[:series_plotindex],
     #         color = py_fillcolor(d),
     #         edgecolor = py_linecolor(d),
     #         linewidth = d[:linewidth],
@@ -594,7 +594,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
     #     end
     #     handle = ax[:hist2d](x, y;
     #         label = d[:label],
-    #         zorder = plt.n,
+    #         zorder = d[:series_plotindex],
     #         bins = d[:bins],
     #         normed = d[:normalize],
     #         weights = d[:weights],
@@ -619,7 +619,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
         end
         handle = ax[:hexbin](x, y;
             label = d[:label],
-            zorder = plt.n,
+            zorder = d[:series_plotindex],
             gridsize = d[:bins],
             linewidths = py_dpi_scale(plt, d[:linewidth]),
             edgecolors = py_linecolor(d),
@@ -660,7 +660,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
         # contour lines
         handle = ax[:contour](x, y, z, levelargs...;
             label = d[:label],
-            zorder = plt.n,
+            zorder = d[:series_plotindex],
             linewidths = py_dpi_scale(plt, d[:linewidth]),
             linestyles = py_linestyle(st, d[:linestyle]),
             cmap = py_linecolormap(d),
@@ -672,7 +672,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
         if d[:fillrange] != nothing
             handle = ax[:contourf](x, y, z, levelargs...;
                 label = d[:label],
-                zorder = plt.n + 0.5,
+                zorder = d[:series_plotindex] + 0.5,
                 cmap = py_fillcolormap(d),
                 extrakw...
             )
@@ -705,7 +705,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
             end
             handle = ax[st == :surface ? :plot_surface : :plot_wireframe](x, y, z;
                 label = d[:label],
-                zorder = plt.n,
+                zorder = d[:series_plotindex],
                 rstride = 1,
                 cstride = 1,
                 linewidth = py_dpi_scale(plt, d[:linewidth]),
@@ -742,7 +742,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
             end
             handle = ax[:plot_trisurf](x, y, z;
                 label = d[:label],
-                zorder = plt.n,
+                zorder = d[:series_plotindex],
                 cmap = py_fillcolormap(d),
                 linewidth = py_dpi_scale(plt, d[:linewidth]),
                 edgecolor = py_linecolor(d),
@@ -766,7 +766,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
             z  # hopefully it's in a data format that will "just work" with imshow
         end
         handle = ax[:imshow](z;
-            zorder = plt.n,
+            zorder = d[:series_plotindex],
             cmap = py_colormap([:black, :white]),
             vmin = 0.0,
             vmax = 1.0
@@ -797,7 +797,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
         
         handle = ax[:pcolormesh](x, y, z;
             label = d[:label],
-            zorder = plt.n,
+            zorder = d[:series_plotindex],
             cmap = py_fillcolormap(d),
             edgecolors = (d[:linewidth] > 0 ? py_linecolor(d) : "face"),
             extrakw...
@@ -820,7 +820,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
         path = py_path(x, y)
         patches = pypatches.pymember("PathPatch")(path;
             label = d[:label],
-            zorder = plt.n,
+            zorder = d[:series_plotindex],
             edgecolor = py_linecolor(d),
             facecolor = py_fillcolor(d),
             linewidth = py_dpi_scale(plt, d[:linewidth]),
@@ -887,7 +887,7 @@ function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
         end
 
         handle = ax[f](args...;
-            zorder = plt.n,
+            zorder = d[:series_plotindex],
             facecolor = py_fillcolor(d),
             linewidths = 0
         )
diff --git a/src/plot.jl b/src/plot.jl
index 663f3003..e10683b7 100644
--- a/src/plot.jl
+++ b/src/plot.jl
@@ -588,75 +588,15 @@ function prepared_object(plt::Plot)
 end
 
 # --------------------------------------------------------------------
+# plot to a Subplot
 
-# function get_indices(orig, labels)
-#     Int[findnext(labels, l, 1) for l in orig]
-# end
-
-# # TODO: remove?? this is the old way of handling discrete data... should be
-# # replaced by the Axis type and logic
-# function setTicksFromStringVector(plt::Plot, d::KW, di::KW, letter)
-#     sym = Symbol(letter)
-#     ticksym = Symbol(letter * "ticks")
-#     pargs = plt.attr
-#     v = di[sym]
-#
-#     # do we really want to do this?
-#     typeof(v) <: AbstractArray || return
-#     isempty(v) && return
-#     trueOrAllTrue(_ -> typeof(_) <: AbstractString, v) || return
-#
-#     # compute the ticks and labels
-#     ticks, labels = if ticksType(pargs[ticksym]) == :ticks_and_labels
-#         # extend the existing ticks and labels. only add to labels if they're new!
-#         ticks, labels = pargs[ticksym]
-#         newlabels = filter(_ -> !(_ in labels), unique(v))
-#         newticks = if isempty(ticks)
-#             collect(1:length(newlabels))
-#         else
-#             maximum(ticks) + collect(1:length(newlabels))
-#         end
-#         ticks = vcat(ticks, newticks)
-#         labels = vcat(labels, newlabels)
-#         ticks, labels
-#     else
-#         # create new ticks and labels
-#         newlabels = unique(v)
-#         collect(1:length(newlabels)), newlabels
-#     end
-#
-#     d[ticksym] = ticks, labels
-#     plt.attr[ticksym] = ticks, labels
-#
-#     # add an origsym field so that later on we can re-compute the x vector if ticks change
-#     origsym = Symbol(letter * "orig")
-#     di[origsym] = v
-#     di[sym] = get_indices(v, labels)
-#
-#     # loop through existing plt.seriesargs and adjust indices if there is an origsym key
-#     for sargs in plt.seriesargs
-#         if haskey(sargs, origsym)
-#             # TODO: might need to call the setindex function instead to trigger a plot update for some backends??
-#             sargs[sym] = get_indices(sargs[origsym], labels)
-#         end
-#     end
-# end
-
-
-# --------------------------------------------------------------------
-
-
-
-# --------------------------------------------------------------------
-
-# function Base.copy(plt::Plot)
-#     backend(plt.backend)
-#     plt2 = plot(; plt.attr...)
-#     for sargs in plt.seriesargs
-#         sargs = filter((k,v) -> haskey(_series_defaults,k), sargs)
-#         plot!(plt2; sargs...)
-#     end
-#     plt2
-# end
+function plot(sp::Subplot, args...; kw...)
+    plt = sp.plt
+    plot(plt, args...; kw..., subplot = findfirst(plt.subplots, sp))
+end
+function plot!(sp::Subplot, args...; kw...)
+    plt = sp.plt
+    plot!(plt, args...; kw..., subplot = findfirst(plt.subplots, sp))
+end
 
 # --------------------------------------------------------------------
diff --git a/src/recipes.jl b/src/recipes.jl
index 3944d8a4..00f55f96 100644
--- a/src/recipes.jl
+++ b/src/recipes.jl
@@ -747,7 +747,14 @@ end
 end
 @deps density path
 
+# ---------------------------------------------------------------------------
+# contourf - filled contours
 
+@recipe function f(::Type{Val{:contourf}}, x, y, z)
+    fillrange := true
+    seriestype := :contour
+    ()
+end
 
 # ---------------------------------------------------------------------------
 # Error Bars
diff --git a/src/series_args.jl b/src/series_args.jl
index ec422e2b..7d491290 100644
--- a/src/series_args.jl
+++ b/src/series_args.jl
@@ -7,7 +7,7 @@
 
 typealias FuncOrFuncs @compat(Union{Function, AVec{Function}})
 
-all3D(d::KW) = trueOrAllTrue(st -> st in (:contour, :heatmap, :surface, :wireframe, :contour3d, :image), get(d, :seriestype, :none))
+all3D(d::KW) = trueOrAllTrue(st -> st in (:contour, :contourf, :heatmap, :surface, :wireframe, :contour3d, :image), get(d, :seriestype, :none))
 
 # missing
 convertToAnyVector(v::@compat(Void), d::KW) = Any[nothing], nothing
diff --git a/src/subplots.jl b/src/subplots.jl
index 06bd3779..6e41bbf5 100644
--- a/src/subplots.jl
+++ b/src/subplots.jl
@@ -39,7 +39,7 @@ function should_add_to_legend(series::Series)
     series.d[:primary] && series.d[:label] != "" &&
         !(series.d[:seriestype] in (
             :hexbin,:histogram2d,:hline,:vline,
-            :contour,:contour3d,:surface,:wireframe,
+            :contour,:contourf,:contour3d,:surface,:wireframe,
             :heatmap, :pie, :image
         ))
 end
diff --git a/src/types.jl b/src/types.jl
index 59bcb5f2..5b97e0e8 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -83,12 +83,13 @@ function Plot()
          Subplot[], false)
 end
 
-# TODO: make a decision... should plt[1] return the first subplot or the first series??
-# Base.getindex(plt::Plot, i::Integer) = plt.subplots[i]
-Base.getindex(plt::Plot, s::Symbol) = plt.spmap[s]
-Base.getindex(plt::Plot, r::Integer, c::Integer) = plt.layout[r,c]
-attr(plt::Plot, k::Symbol) = plt.attr[k]
-attr!(plt::Plot, v, k::Symbol) = (plt.attr[k] = v)
+# -----------------------------------------------------------------------
 
+Base.getindex(plt::Plot, i::Integer) = plt.subplots[i]
+Base.getindex(plt::Plot, r::Integer, c::Integer) = plt.layout[r,c]
+# attr(plt::Plot, k::Symbol) = plt.attr[k]
+# attr!(plt::Plot, v, k::Symbol) = (plt.attr[k] = v)
+
+Base.getindex(sp::Subplot, i::Integer) = series_list(sp)[i]
 
 # -----------------------------------------------------------------------
diff --git a/src/utils.jl b/src/utils.jl
index 67fea26e..b8add3d0 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -642,7 +642,7 @@ end
 # -------------------------------------------------------
 # indexing notation
 
-Base.getindex(plt::Plot, i::Integer) = getxy(plt, i)
+# Base.getindex(plt::Plot, i::Integer) = getxy(plt, i)
 Base.setindex!{X,Y}(plt::Plot, xy::Tuple{X,Y}, i::Integer) = setxy!(plt, xy, i)
 Base.setindex!{X,Y,Z}(plt::Plot, xyz::Tuple{X,Y,Z}, i::Integer) = setxyz!(plt, xyz, i)