diff --git a/deps/agg/include/agg_color_gray.h b/deps/agg/include/agg_color_gray.h
index 19d31d38e..2fc01bd71 100644
--- a/deps/agg/include/agg_color_gray.h
+++ b/deps/agg/include/agg_color_gray.h
@@ -746,6 +746,13 @@ struct gray32
     value_type v;
     value_type a;
 
+    enum base_scale_e
+    {
+        base_shift = 8,
+        base_scale = 1 << base_shift,
+        base_mask  = base_scale - 1,
+    };
+
     // Calculate grayscale value as per ITU-R BT.709.
     static value_type luminance(double r, double g, double b)
     {
diff --git a/deps/agg/include/agg_pixfmt_base.h b/deps/agg/include/agg_pixfmt_base.h
new file mode 100644
index 000000000..57ae19cfe
--- /dev/null
+++ b/deps/agg/include/agg_pixfmt_base.h
@@ -0,0 +1,97 @@
+//----------------------------------------------------------------------------
+// Anti-Grain Geometry - Version 2.4
+// Copyright (C) 2002-2005 Maxim Shemanarev (http://www.antigrain.com)
+//
+// Permission to copy, use, modify, sell and distribute this software 
+// is granted provided this copyright notice appears in all copies. 
+// This software is provided "as is" without express or implied
+// warranty, and with no claim as to its suitability for any purpose.
+//
+//----------------------------------------------------------------------------
+// Contact: mcseem@antigrain.com
+//          mcseemagg@yahoo.com
+//          http://www.antigrain.com
+//----------------------------------------------------------------------------
+
+#ifndef AGG_PIXFMT_BASE_INCLUDED
+#define AGG_PIXFMT_BASE_INCLUDED
+
+#include "agg_basics.h"
+#include "agg_color_gray.h"
+#include "agg_color_rgba.h"
+
+namespace agg
+{
+    struct pixfmt_gray_tag
+    {
+    };
+
+    struct pixfmt_rgb_tag
+    {
+    };
+
+    struct pixfmt_rgba_tag
+    {
+    };
+
+    //--------------------------------------------------------------blender_base
+    template<class ColorT, class Order = void> 
+    struct blender_base
+    {
+        typedef ColorT color_type;
+        typedef Order order_type;
+        typedef typename color_type::value_type value_type;
+
+        static rgba get(value_type r, value_type g, value_type b, value_type a, cover_type cover = cover_full)
+        {
+            if (cover > cover_none)
+            {
+                rgba c(
+                    color_type::to_double(r), 
+                    color_type::to_double(g), 
+                    color_type::to_double(b), 
+                    color_type::to_double(a));
+
+                if (cover < cover_full)
+                {
+                    double x = double(cover) / cover_full;
+                    c.r *= x;
+                    c.g *= x;
+                    c.b *= x;
+                    c.a *= x;
+                }
+
+                return c;
+            }
+            else return rgba::no_color();
+        }
+
+        static rgba get(const value_type* p, cover_type cover = cover_full)
+        {
+            return get(
+                p[order_type::R], 
+                p[order_type::G], 
+                p[order_type::B], 
+                p[order_type::A], 
+                cover);
+        }
+
+        static void set(value_type* p, value_type r, value_type g, value_type b, value_type a)
+        {
+            p[order_type::R] = r;
+            p[order_type::G] = g;
+            p[order_type::B] = b;
+            p[order_type::A] = a;
+        }
+
+        static void set(value_type* p, const rgba& c)
+        {
+            p[order_type::R] = color_type::from_double(c.r);
+            p[order_type::G] = color_type::from_double(c.g);
+            p[order_type::B] = color_type::from_double(c.b);
+            p[order_type::A] = color_type::from_double(c.a);
+        }
+    };
+}
+
+#endif
diff --git a/include/mapnik/image_compositing.hpp b/include/mapnik/image_compositing.hpp
index aab3fcc2c..a769442ac 100644
--- a/include/mapnik/image_compositing.hpp
+++ b/include/mapnik/image_compositing.hpp
@@ -90,13 +90,21 @@ MAPNIK_DECL void composite(T1 & dst, T2 & src,
                            int dy=0,
                            bool premultiply_src=false);
 
-extern template MAPNIK_DECL void composite<mapnik::image_data_32,mapnik::image_data_32>(mapnik::image_data_32 & dst,
-                           mapnik::image_data_32 & src,
-                           composite_mode_e mode,
-                           float opacity,
-                           int dx,
-                           int dy,
-                           bool premultiply_src);
+//extern template MAPNIK_DECL void composite<mapnik::image_data_32,mapnik::image_data_32>(mapnik::image_data_32 & dst,
+//                           mapnik::image_data_32 & src,
+//                           composite_mode_e mode,
+//                           float opacity,
+//                           int dx,
+//                           int dy,
+//                           bool premultiply_src);
+
+//extern template MAPNIK_DECL void composite<mapnik::image_data_float32,mapnik::image_data_float32>(mapnik::image_data_float32 & dst,
+//                           mapnik::image_data_float32 & src,
+//                           composite_mode_e mode,
+//                           float opacity,
+//                           int dx,
+//                           int dy,
+//                           bool premultiply_src);
 
 }
 #endif // MAPNIK_IMAGE_COMPOSITING_HPP
diff --git a/include/mapnik/image_data.hpp b/include/mapnik/image_data.hpp
index 78b4acccb..64dbbb66e 100644
--- a/include/mapnik/image_data.hpp
+++ b/include/mapnik/image_data.hpp
@@ -190,7 +190,8 @@ private:
 };
 
 using image_data_32 = image_data<std::uint32_t>;
-using image_data_8 = image_data<byte> ;
+using image_data_8 = image_data<std::uint8_t> ;
+using image_data_16 = image_data<std::uint16_t>;
 using image_data_float32 = image_data<float>;
 }
 
diff --git a/include/mapnik/raster.hpp b/include/mapnik/raster.hpp
index 791052de5..a0c7de07d 100644
--- a/include/mapnik/raster.hpp
+++ b/include/mapnik/raster.hpp
@@ -33,7 +33,7 @@
 
 namespace mapnik {
 
-using image_data_base = util::variant<image_data_32, image_data_8, image_data_float32>;
+using image_data_base = util::variant<image_data_32, image_data_8, image_data_16, image_data_float32>;
 
 namespace detail {
 
diff --git a/include/mapnik/raster_colorizer.hpp b/include/mapnik/raster_colorizer.hpp
index 358383b7b..530d1e703 100644
--- a/include/mapnik/raster_colorizer.hpp
+++ b/include/mapnik/raster_colorizer.hpp
@@ -199,7 +199,7 @@ public:
     //! \brief Colorize a raster
     //!
     //! \param[in, out] raster A raster stored in float32 single channel format, which gets colorized in place.
-    void colorize(std::shared_ptr<raster> const& raster, feature_impl const& f) const;
+    void colorize(raster & ras, feature_impl const& f) const;
 
 
     //! \brief Perform the translation of input to output
diff --git a/include/mapnik/renderer_common/process_raster_symbolizer.hpp b/include/mapnik/renderer_common/process_raster_symbolizer.hpp
index dcc956e0d..56b75661f 100644
--- a/include/mapnik/renderer_common/process_raster_symbolizer.hpp
+++ b/include/mapnik/renderer_common/process_raster_symbolizer.hpp
@@ -34,6 +34,10 @@
 // agg
 #include "agg_rendering_buffer.h"
 #include "agg_pixfmt_rgba.h"
+#include "agg_pixfmt_gray.h"
+#include "agg_rasterizer_scanline_aa.h"
+#include "agg_scanline_u.h"
+#include "agg_renderer_scanline.h"
 
 namespace mapnik {
 
@@ -47,11 +51,6 @@ void render_raster_symbolizer(raster_symbolizer const &sym,
     raster_ptr const& source = feature.get_raster();
     if (source)
     {
-        // If there's a colorizer defined, use it to color the raster in-place
-        raster_colorizer_ptr colorizer = get<raster_colorizer_ptr>(sym, keys::colorizer);
-        if (colorizer)
-            colorizer->colorize(source,feature);
-
         box2d<double> target_ext = box2d<double>(source->ext_);
         prj_trans.backward(target_ext, PROJ_ENVELOPE_POINTS);
         box2d<double> ext = common.t_.forward(target_ext);
@@ -129,8 +128,7 @@ void render_raster_symbolizer(raster_symbolizer const &sym,
                     if (source->data_.is<image_data_32>())
                     {
                         image_data_32 data(raster_width, raster_height);
-                        raster target(target_ext, data, source->get_filter_factor());
-                        scale_image_agg<image_data_32>(util::get<image_data_32>(target.data_),
+                        scale_image_agg<image_data_32>(data,
                                                        util::get<image_data_32>(source->data_),
                                                        scaling_method,
                                                        image_ratio_x,
@@ -138,23 +136,37 @@ void render_raster_symbolizer(raster_symbolizer const &sym,
                                                        0.0,
                                                        0.0,
                                                        source->get_filter_factor());
-                        composite(util::get<image_data_32>(target.data_), comp_op, opacity, start_x, start_y);
+                        composite(data, comp_op, opacity, start_x, start_y);
                     }
                     else if (source->data_.is<image_data_float32>())
                     {
                         std::cerr << "#3 source->data float32" << std::endl;
-                        image_data_float32 data(raster_width, raster_height);
-                        raster target(target_ext, data, source->get_filter_factor());
-                        //scale_image_agg<image_data_float32>(util::get<image_data_float32>(target.data_),
-                        //                                    util::get<image_data_float32>(source->data_),
-                        //                                    scaling_method,
-                        //                                    image_ratio_x,
-                        //                                    image_ratio_y,
-                        //                                    0.0,
-                        //                                    0.0,
-                        //                                    source->get_filter_factor());
-                        // composite is no-op
+                        raster target(target_ext, image_data_float32(raster_width, raster_height), source->get_filter_factor());
+                        scale_image_agg<image_data_float32>(util::get<image_data_float32>(target.data_),
+                                                            util::get<image_data_float32>(source->data_),
+                                                            scaling_method,
+                                                            image_ratio_x,
+                                                            image_ratio_y,
+                                                            0.0,
+                                                            0.0,
+                                                            source->get_filter_factor());
 
+                        image_data_float32 & data = util::get<image_data_float32>(target.data_);
+                        image_data_32 dst(raster_width, raster_height);
+                        for (unsigned x = 0; x < dst.width(); ++x)
+                        {
+                            for (unsigned y = 0; y < dst.height(); ++y)
+                            {
+                                //std::cerr << data(x,y) << std::endl;
+                                unsigned val = unsigned(255*data(x,y)+0.5);
+                                dst(x,y) = color(val,val,val, unsigned(opacity*255+0.5)).rgba();
+                            }
+                        }
+                        raster out(target_ext, dst,source->get_filter_factor());
+                        // If there's a colorizer defined, use it to color the raster in-place
+                        raster_colorizer_ptr colorizer = get<raster_colorizer_ptr>(sym, keys::colorizer);
+                        if (colorizer) colorizer->colorize(out,feature);
+                        composite(util::get<image_data_32>(out.data_), comp_op, opacity, start_x, start_y);
                     }
                 }
             }
diff --git a/src/image_compositing.cpp b/src/image_compositing.cpp
index 1bca16274..73d64cb0e 100644
--- a/src/image_compositing.cpp
+++ b/src/image_compositing.cpp
@@ -38,8 +38,10 @@
 #include "agg_scanline_u.h"
 #include "agg_renderer_scanline.h"
 #include "agg_pixfmt_rgba.h"
+#include "agg_pixfmt_gray.h"
 #include "agg_color_rgba.h"
 
+
 namespace mapnik
 {
 
@@ -121,8 +123,8 @@ For example, if you generate some pattern with AGG (premultiplied) and would lik
 */
 
 
-template <typename T1, typename T2>
-void composite(T1 & dst, T2 & src, composite_mode_e mode,
+template <>
+void composite(mapnik::image_data_32 & dst, mapnik::image_data_32 & src, composite_mode_e mode,
                float opacity,
                int dx,
                int dy,
@@ -146,6 +148,32 @@ void composite(T1 & dst, T2 & src, composite_mode_e mode,
     ren.blend_from(pixf_mask,0,dx,dy,unsigned(255*opacity));
 }
 
+template <>
+void composite(mapnik::image_data_float32 & dst, mapnik::image_data_float32 & src, composite_mode_e mode,
+               float opacity,
+               int dx,
+               int dy,
+               bool premultiply_src)
+{
+    //using color = agg::gray32;
+    //using order = agg::order_gray32;
+    //using blender_type = agg::comp_op_adaptor_rgba_pre<color, order>;
+    using pixfmt_type = agg::pixfmt_gray32;//agg::pixfmt_custom_blend_rgba<blender_type, agg::rendering_buffer>;
+    using renderer_type = agg::renderer_base<pixfmt_type>;
+
+    agg::rendering_buffer dst_buffer(dst.getBytes(),dst.width(),dst.height(),dst.width());
+    agg::rendering_buffer src_buffer(src.getBytes(),src.width(),src.height(),src.width());
+
+    pixfmt_type pixf(dst_buffer);
+    //pixf.comp_op(static_cast<agg::comp_op_e>(mode));
+
+    agg::pixfmt_gray32 pixf_mask(src_buffer);
+    //if (premultiply_src)  pixf_mask.premultiply();
+    renderer_type ren(pixf);
+    //ren.blend_from(pixf_mask,0,dx,dy,agg::cover_full);//unsigned(255*opacity));
+    ren.copy_from(pixf_mask,0,dx,dy);//unsigned(255*opacity));
+}
+
 template void composite<mapnik::image_data_32,mapnik::image_data_32>(mapnik::image_data_32&,
                                                                      mapnik::image_data_32&,
                                                                      composite_mode_e,
diff --git a/src/image_scaling.cpp b/src/image_scaling.cpp
index 43c2b8a8d..545910f03 100644
--- a/src/image_scaling.cpp
+++ b/src/image_scaling.cpp
@@ -37,12 +37,14 @@
 // agg
 #include "agg_image_accessors.h"
 #include "agg_pixfmt_rgba.h"
+#include "agg_pixfmt_gray.h"
 #include "agg_color_rgba.h"
 #include "agg_rasterizer_scanline_aa.h"
 #include "agg_renderer_scanline.h"
 #include "agg_rendering_buffer.h"
 #include "agg_scanline_u.h"
 #include "agg_span_allocator.h"
+#include "agg_span_image_filter_gray.h"
 #include "agg_span_image_filter_rgba.h"
 #include "agg_span_interpolator_linear.h"
 #include "agg_trans_affine.h"
@@ -94,14 +96,15 @@ boost::optional<std::string> scaling_method_to_string(scaling_method_e scaling_m
     return mode;
 }
 
-template <> void scale_image_agg<image_data_32>(image_data_32 & target,
-                                                image_data_32 const& source,
-                                                scaling_method_e scaling_method,
-                                                double image_ratio_x,
-                                                double image_ratio_y,
-                                                double x_off_f,
-                                                double y_off_f,
-                                                double filter_factor)
+template <>
+void scale_image_agg<image_data_32>(image_data_32 & target,
+                                    image_data_32 const& source,
+                                    scaling_method_e scaling_method,
+                                    double image_ratio_x,
+                                    double image_ratio_y,
+                                    double x_off_f,
+                                    double y_off_f,
+                                    double filter_factor)
 {
     // "the image filters should work namely in the premultiplied color space"
     // http://old.nabble.com/Re:--AGG--Basic-image-transformations-p1110665.html
@@ -204,13 +207,108 @@ template <> void scale_image_agg<image_data_32>(image_data_32 & target,
     agg::render_scanlines_aa(ras, sl, rb_dst_pre, sa, sg);
 }
 
-template void scale_image_agg<image_data_32>(image_data_32& target,
-                                             const image_data_32& source,
-                                             scaling_method_e scaling_method,
-                                             double image_ratio_x,
-                                             double image_ratio_y,
-                                             double x_off_f,
-                                             double y_off_f,
-                                             double filter_factor);
+
+template <>
+void scale_image_agg<image_data_float32>(image_data_float32 & target,
+                                    image_data_float32 const& source,
+                                    scaling_method_e scaling_method,
+                                    double image_ratio_x,
+                                    double image_ratio_y,
+                                    double x_off_f,
+                                    double y_off_f,
+                                    double filter_factor)
+{
+    using pixfmt_pre = agg::pixfmt_gray32_pre;
+    using renderer_base_pre = agg::renderer_base<pixfmt_pre>;
+
+    // define some stuff we'll use soon
+    agg::rasterizer_scanline_aa<> ras;
+    agg::scanline_u8 sl;
+    agg::span_allocator<agg::gray32> sa;
+    agg::image_filter_lut filter;
+
+    // initialize source AGG buffer
+    agg::rendering_buffer rbuf_src(const_cast<unsigned char*>(source.getBytes()), source.width(), source.height(), source.width() * 4);
+    pixfmt_pre pixf_src(rbuf_src);
+    using img_src_type = agg::image_accessor_clone<pixfmt_pre>;
+    img_src_type img_src(pixf_src);
+
+    // initialize destination AGG buffer (with transparency)
+    agg::rendering_buffer rbuf_dst(target.getBytes(), target.width(), target.height(), target.width() * 4);
+    pixfmt_pre pixf_dst(rbuf_dst);
+    renderer_base_pre rb_dst_pre(pixf_dst);
+
+    // create a scaling matrix
+    agg::trans_affine img_mtx;
+    img_mtx /= agg::trans_affine_scaling(image_ratio_x, image_ratio_y);
+
+    // create a linear interpolator for our scaling matrix
+    using interpolator_type = agg::span_interpolator_linear<>;
+    interpolator_type interpolator(img_mtx);
+
+    // draw an anticlockwise polygon to render our image into
+    double scaled_width = target.width();
+    double scaled_height = target.height();
+    ras.reset();
+    ras.move_to_d(x_off_f,                y_off_f);
+    ras.line_to_d(x_off_f + scaled_width, y_off_f);
+    ras.line_to_d(x_off_f + scaled_width, y_off_f + scaled_height);
+    ras.line_to_d(x_off_f,                y_off_f + scaled_height);
+
+    switch(scaling_method)
+    {
+    case SCALING_NEAR:
+    {
+        using span_gen_type = agg::span_image_filter_gray_nn<img_src_type, interpolator_type>;
+        span_gen_type sg(img_src, interpolator);
+        agg::render_scanlines_aa(ras, sl, rb_dst_pre, sa, sg);
+        return;
+    }
+    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(filter_factor), true); break;
+    case SCALING_LANCZOS:
+        filter.calculate(agg::image_filter_lanczos(filter_factor), true); break;
+    case SCALING_BLACKMAN:
+        filter.calculate(agg::image_filter_blackman(filter_factor), true); break;
+    }
+    using span_gen_type = agg::span_image_resample_gray_affine<img_src_type>;
+    span_gen_type sg(img_src, interpolator, filter);
+    agg::render_scanlines_aa(ras, sl, rb_dst_pre, sa, sg);
+}
+
+//template void scale_image_agg<image_data_32>(image_data_32& target,
+//                                             const image_data_32& source,
+//                                             scaling_method_e scaling_method,
+//                                             double image_ratio_x,
+//                                             double image_ratio_y,
+//                                             double x_off_f,
+//                                             double y_off_f,
+//                                             double filter_factor);
 
 }
diff --git a/src/raster_colorizer.cpp b/src/raster_colorizer.cpp
index 3cf0b0850..8fd8634f6 100644
--- a/src/raster_colorizer.cpp
+++ b/src/raster_colorizer.cpp
@@ -122,16 +122,16 @@ bool raster_colorizer::add_stop(colorizer_stop const& stop)
     return true;
 }
 
-void raster_colorizer::colorize(raster_ptr const& raster, feature_impl const& f) const
+void raster_colorizer::colorize(raster & ras, feature_impl const& f) const
 {
-    unsigned* imageData = reinterpret_cast<unsigned*>(raster->data_.getBytes());
+    unsigned * imageData = reinterpret_cast<unsigned*>(ras.data_.getBytes());
 
-    int len = raster->data_.width() * raster->data_.height();
-    boost::optional<double> const& nodata = raster->nodata();
+    int len = ras.data_.width() * ras.data_.height();
+    boost::optional<double> const& nodata = ras.nodata();
     for (int i=0; i<len; ++i)
     {
         // the GDAL plugin reads single bands as floats
-        float value = *reinterpret_cast<float*> (&imageData[i]);
+        float value = *reinterpret_cast<float *> (&imageData[i]);
         if (nodata && (std::fabs(value - *nodata) < epsilon_))
         {
             imageData[i] = 0;
@@ -153,7 +153,7 @@ unsigned raster_colorizer::get_color(float value) const
     int stopCount = stops_.size();
 
     //use default color if no stops
-    if(stopCount == 0)
+    if (stopCount == 0)
     {
         return default_color_.rgba();
     }