diff --git a/include/mapnik/warp.hpp b/include/mapnik/warp.hpp
index 0c3f13608..bd43307f1 100644
--- a/include/mapnik/warp.hpp
+++ b/include/mapnik/warp.hpp
@@ -33,11 +33,11 @@ class raster;
 class proj_transform;
 
 MAPNIK_DECL void reproject_and_scale_raster(raster & target,
-                                raster const& source,
-                                proj_transform const& prj_trans,
-                                double offset_x, double offset_y,
-                                unsigned mesh_size,
-                                scaling_method_e scaling_method);
+                                            raster const& source,
+                                            proj_transform const& prj_trans,
+                                            double offset_x, double offset_y,
+                                            unsigned mesh_size,
+                                            scaling_method_e scaling_method);
 
 }
 
diff --git a/src/warp.cpp b/src/warp.cpp
index eaf93c93d..b66a5180f 100644
--- a/src/warp.cpp
+++ b/src/warp.cpp
@@ -24,6 +24,7 @@
 #include <mapnik/warp.hpp>
 #include <mapnik/config.hpp>
 #include <mapnik/image_data.hpp>
+#include <mapnik/image_scaling_traits.hpp>
 #include <mapnik/image_util.hpp>
 #include <mapnik/box2d.hpp>
 #include <mapnik/view_transform.hpp>
@@ -47,19 +48,29 @@
 
 namespace mapnik {
 
-void reproject_and_scale_raster(raster & target, raster const& source,
-                                proj_transform const& prj_trans,
-                                double offset_x, double offset_y,
-                                unsigned mesh_size,
-                                scaling_method_e scaling_method)
-{
-    view_transform ts(source.data_.width(), source.data_.height(),
-                      source.ext_);
-    view_transform tt(target.data_.width(), target.data_.height(),
-                      target.ext_, offset_x, offset_y);
+namespace  detail {
 
-    std::size_t mesh_nx = std::ceil(source.data_.width()/double(mesh_size) + 1);
-    std::size_t mesh_ny = std::ceil(source.data_.height()/double(mesh_size) + 1);
+template <typename T>
+void warp_image (T & target, T const& source, proj_transform const& prj_trans,
+                 box2d<double> const& target_ext, box2d<double> const& source_ext,
+                 double offset_x, double offset_y, unsigned mesh_size, scaling_method_e scaling_method, double filter_factor)
+{
+    using image_data_type = T;
+    using pixel_type = typename image_data_type::pixel_type;
+    using pixfmt_pre = typename detail::agg_scaling_traits<image_data_type>::pixfmt_pre;
+    using color_type = typename detail::agg_scaling_traits<image_data_type>::color_type;
+    using renderer_base = agg::renderer_base<pixfmt_pre>;
+    using interpolator_type = typename detail::agg_scaling_traits<image_data_type>::interpolator_type;
+
+    constexpr std::size_t pixel_size = sizeof(pixel_type);
+
+    view_transform ts(source.width(), source.height(),
+                      source_ext);
+    view_transform tt(target.width(), target.height(),
+                      target_ext, offset_x, offset_y);
+
+    std::size_t mesh_nx = std::ceil(source.width()/double(mesh_size) + 1);
+    std::size_t mesh_ny = std::ceil(source.height()/double(mesh_size) + 1);
 
     image_data<double> xs(mesh_nx, mesh_ny);
     image_data<double> ys(mesh_nx, mesh_ny);
@@ -69,81 +80,36 @@ void reproject_and_scale_raster(raster & target, raster const& source,
     {
         for (std::size_t i=0; i<mesh_nx; ++i)
         {
-            xs(i,j) = std::min(i*mesh_size,source.data_.width());
-            ys(i,j) = std::min(j*mesh_size,source.data_.height());
+            xs(i,j) = std::min(i*mesh_size,source.width());
+            ys(i,j) = std::min(j*mesh_size,source.height());
             ts.backward(&xs(i,j), &ys(i,j));
         }
     }
     prj_trans.backward(xs.getData(), ys.getData(), nullptr, mesh_nx*mesh_ny);
 
-    // Initialize AGG objects
-    using pixfmt = agg::pixfmt_rgba32_pre;
-    using color_type = pixfmt::color_type;
-    using renderer_base = agg::renderer_base<pixfmt>;
-
     agg::rasterizer_scanline_aa<> rasterizer;
     agg::scanline_bin scanline;
-    agg::rendering_buffer buf(target.data_.getBytes(),
-                              target.data_.width(),
-                              target.data_.height(),
-                              target.data_.width()*4);
-    pixfmt pixf(buf);
+    agg::rendering_buffer buf(target.getBytes(),
+                              target.width(),
+                              target.height(),
+                              target.width() * pixel_size);
+    pixfmt_pre pixf(buf);
     renderer_base rb(pixf);
-    rasterizer.clip_box(0, 0, target.data_.width(), target.data_.height());
+    rasterizer.clip_box(0, 0, target.width(), target.height());
     agg::rendering_buffer buf_tile(
-        const_cast<unsigned char*>(source.data_.getBytes()),
-        source.data_.width(),
-        source.data_.height(),
-        source.data_.width() * 4);
+        const_cast<unsigned char*>(source.getBytes()),
+        source.width(),
+        source.height(),
+        source.width() * pixel_size);
 
-    pixfmt pixf_tile(buf_tile);
+    pixfmt_pre pixf_tile(buf_tile);
 
-    using img_accessor_type = agg::image_accessor_clone<pixfmt>;
+    using img_accessor_type = agg::image_accessor_clone<pixfmt_pre>;
     img_accessor_type ia(pixf_tile);
 
     agg::span_allocator<color_type> sa;
-
-    // Initialize filter
-    agg::image_filter_lut filter;
-    switch(scaling_method)
-    {
-    case SCALING_NEAR: break;
-    case SCALING_BILINEAR:
-        filter.calculate(agg::image_filter_bilinear(), true); break;
-    case SCALING_BICUBIC:
-        filter.calculate(agg::image_filter_bicubic(), true); break;
-    case SCALING_SPLINE16:
-        filter.calculate(agg::image_filter_spline16(), true); break;
-    case SCALING_SPLINE36:
-        filter.calculate(agg::image_filter_spline36(), true); break;
-    case SCALING_HANNING:
-        filter.calculate(agg::image_filter_hanning(), true); break;
-    case SCALING_HAMMING:
-        filter.calculate(agg::image_filter_hamming(), true); break;
-    case SCALING_HERMITE:
-        filter.calculate(agg::image_filter_hermite(), true); break;
-    case SCALING_KAISER:
-        filter.calculate(agg::image_filter_kaiser(), true); break;
-    case SCALING_QUADRIC:
-        filter.calculate(agg::image_filter_quadric(), true); break;
-    case SCALING_CATROM:
-        filter.calculate(agg::image_filter_catrom(), true); break;
-    case SCALING_GAUSSIAN:
-        filter.calculate(agg::image_filter_gaussian(), true); break;
-    case SCALING_BESSEL:
-        filter.calculate(agg::image_filter_bessel(), true); break;
-    case SCALING_MITCHELL:
-        filter.calculate(agg::image_filter_mitchell(), true); break;
-    case SCALING_SINC:
-        filter.calculate(agg::image_filter_sinc(source.get_filter_factor()), true); break;
-    case SCALING_LANCZOS:
-        filter.calculate(agg::image_filter_lanczos(source.get_filter_factor()), true); break;
-    case SCALING_BLACKMAN:
-        filter.calculate(agg::image_filter_blackman(source.get_filter_factor()), true); break;
-    }
-
     // Project mesh cells into target interpolating raster inside each one
-    for(std::size_t j = 0; j < mesh_ny - 1; ++j)
+    for (std::size_t j = 0; j < mesh_ny - 1; ++j)
     {
         for (std::size_t i = 0; i < mesh_nx - 1; ++i)
         {
@@ -166,33 +132,82 @@ void reproject_and_scale_raster(raster & target, raster const& source,
             std::size_t y0 = j * mesh_size;
             std::size_t x1 = (i+1) * mesh_size;
             std::size_t y1 = (j+1) * mesh_size;
-            x1 = std::min(x1, source.data_.width());
-            y1 = std::min(y1, source.data_.height());
+            x1 = std::min(x1, source.width());
+            y1 = std::min(y1, source.height());
             agg::trans_affine tr(polygon, x0, y0, x1, y1);
             if (tr.is_valid())
             {
-                using interpolator_type = agg::span_interpolator_linear<agg::trans_affine>;
                 interpolator_type interpolator(tr);
-
                 if (scaling_method == SCALING_NEAR)
                 {
-                    using span_gen_type = agg::span_image_filter_rgba_nn
-                        <img_accessor_type, interpolator_type>;
+                    using span_gen_type = typename detail::agg_scaling_traits<image_data_type>::span_image_filter;
                     span_gen_type sg(ia, interpolator);
-                    agg::render_scanlines_bin(rasterizer, scanline, rb,
-                                             sa, sg);
+                    agg::render_scanlines_bin(rasterizer, scanline, rb, sa, sg);
                 }
                 else
                 {
-                    using span_gen_type = agg::span_image_resample_rgba_affine
-                        <img_accessor_type>;
+                    using span_gen_type = typename detail::agg_scaling_traits<image_data_type>::span_image_resample_affine;
+                    agg::image_filter_lut filter;
+                    set_scaling_method(filter, scaling_method, filter_factor);
                     span_gen_type sg(ia, interpolator, filter);
-                    agg::render_scanlines_bin(rasterizer, scanline, rb,
-                                             sa, sg);
+                    agg::render_scanlines_bin(rasterizer, scanline, rb, sa, sg);
                 }
             }
 
         }
     }
 }
+
+struct warp_image_visitor : util::static_visitor<void>
+{
+    warp_image_visitor (raster & target_raster, proj_transform const& prj_trans, box2d<double> const& source_ext,
+                        double offset_x, double offset_y, unsigned mesh_size,
+                        scaling_method_e scaling_method, double filter_factor)
+        : target_raster_(target_raster),
+          prj_trans_(prj_trans),
+          source_ext_(source_ext),
+          offset_x_(offset_x),
+          offset_y_(offset_y),
+          mesh_size_(mesh_size),
+          scaling_method_(scaling_method),
+        filter_factor_(filter_factor) {}
+
+    void operator() (image_data_null const&) {}
+
+    template <typename T>
+    void operator() (T const& source)
+    {
+        using image_data_type = T;
+        //source and target image data types must match
+        if (target_raster_.data_.template is<image_data_type>())
+        {
+            image_data_type & target = util::get<image_data_type>(target_raster_.data_);
+            detail::warp_image (target, source, prj_trans_, target_raster_.ext_, source_ext_,
+                                offset_x_, offset_y_, mesh_size_, scaling_method_, filter_factor_);
+        }
+    }
+
+    raster & target_raster_;
+    proj_transform const& prj_trans_;
+    box2d<double> const& source_ext_;
+    double offset_x_;
+    double offset_y_;
+    unsigned mesh_size_;
+    scaling_method_e scaling_method_;
+    double filter_factor_;
+};
+
+}
+
+void reproject_and_scale_raster(raster & target, raster const& source,
+                                proj_transform const& prj_trans,
+                                double offset_x, double offset_y,
+                                unsigned mesh_size,
+                                scaling_method_e scaling_method)
+{
+    detail::warp_image_visitor warper(target, prj_trans, source.ext_, offset_x, offset_y, mesh_size,
+                                      scaling_method, source.get_filter_factor());
+    util::apply_visitor(warper, source.data_);
+}
+
 }// namespace mapnik