diff --git a/include/mapnik/image_util.hpp b/include/mapnik/image_util.hpp
index 22e6574ae..83f09ba23 100644
--- a/include/mapnik/image_util.hpp
+++ b/include/mapnik/image_util.hpp
@@ -218,7 +218,7 @@ inline void scale_image (Image& target,const Image& source)
 }
 
 template <typename Image>
-inline void scale_image_bilinear (Image& target,const Image& source)
+inline void scale_image_bilinear (Image& target,const Image& source, double x_off_f=0, double y_off_f=0)
 {
 
     int source_width=source.width();
@@ -232,28 +232,43 @@ inline void scale_image_bilinear (Image& target,const Image& source)
     int x=0,y=0,xs=0,ys=0;
     int tw2 = target_width/2;
     int th2 = target_height/2;
+    int offs_x = int(round((source_width-target_width-x_off_f*2*source_width)/2));
+    int offs_y = int(round((source_height-target_height-y_off_f*2*source_height)/2));
+    unsigned yprt, yprt1, xprt, xprt1;
 
+    //no scaling or subpixel offset
+    if (target_height == source_height && target_width == source_width && offs_x == 0 && offs_y == 0){
+        for (y=0;y<target_height;++y)
+            target.setRow(y,source.getRow(y),target_width);
+        return;
+    }
 
     for (y=0;y<target_height;++y)
     {
-        ys = y*source_height/target_height;
+        ys = (y*source_height+offs_y)/target_height;
         int ys1 = ys+1;
         if (ys1>=source_height)
             ys1--;
-        unsigned yprt = y*source_height%target_height;
-        unsigned yprt1 = target_height-yprt;
+        if (ys<0)
+            ys=ys1=0;
+        if (source_height/2<target_height)
+            yprt = (y*source_height+offs_y)%target_height;
+        else
+            yprt = th2;
+        yprt1 = target_height-yprt;
         for (x=0;x<target_width;++x)
         {
-            xs = x*source_width/target_width;
-            if (source_width>=target_width || source_height>=target_height){
-                target(x,y)=source(xs,ys);
-                continue;
-            }
-            unsigned xprt = x*source_width%target_width;
-            unsigned xprt1 = target_width-xprt;
+            xs = (x*source_width+offs_x)/target_width;
+            if (source_width/2<target_width)
+                xprt = (x*source_width+offs_x)%target_width;
+            else
+                xprt = tw2;
+            xprt1 = target_width-xprt;
             int xs1 = xs+1;
             if (xs1>=source_width)
                 xs1--;
+            if (xs<0)
+                xs=xs1=0;
 
             unsigned a = source(xs,ys);
             unsigned b = source(xs1,ys);
@@ -290,7 +305,7 @@ inline void scale_image_bilinear (Image& target,const Image& source)
 }
 
 template <typename Image>
-inline void scale_image_bilinear8 (Image& target,const Image& source)
+inline void scale_image_bilinear8 (Image& target,const Image& source, double x_off_f=0, double y_off_f=0)
 {
 
     int source_width=source.width();
@@ -304,28 +319,43 @@ inline void scale_image_bilinear8 (Image& target,const Image& source)
     int x=0,y=0,xs=0,ys=0;
     int tw2 = target_width/2;
     int th2 = target_height/2;
+    int offs_x = int(round((source_width-target_width-x_off_f*2*source_width)/2));
+    int offs_y = int(round((source_height-target_height-y_off_f*2*source_height)/2));
+    unsigned yprt, yprt1, xprt, xprt1;
 
+    //no scaling or subpixel offset
+    if (target_height == source_height && target_width == source_width && offs_x == 0 && offs_y == 0){
+        for (y=0;y<target_height;++y)
+            target.setRow(y,source.getRow(y),target_width);
+        return;
+    }
 
     for (y=0;y<target_height;++y)
     {
-        ys = y*source_height/target_height;
+        ys = (y*source_height+offs_y)/target_height;
         int ys1 = ys+1;
         if (ys1>=source_height)
             ys1--;
-        unsigned yprt = y*source_height%target_height;
-        unsigned yprt1 = target_height-yprt;
+        if (ys<0)
+            ys=ys1=0;
+        if (source_height/2<target_height)
+            yprt = (y*source_height+offs_y)%target_height;
+        else
+            yprt = th2;
+        yprt1 = target_height-yprt;
         for (x=0;x<target_width;++x)
         {
-            xs = x*source_width/target_width;
-            if (source_width>=target_width || source_height>=target_height){
-                target(x,y)=source(xs,ys);
-                continue;
-            }
-            unsigned xprt = x*source_width%target_width;
-            unsigned xprt1 = target_width-xprt;
+            xs = (x*source_width+offs_x)/target_width;
+            if (source_width/2<target_width)
+                xprt = (x*source_width+offs_x)%target_width;
+            else
+                xprt = tw2;
+            xprt1 = target_width-xprt;
             int xs1 = xs+1;
             if (xs1>=source_width)
                 xs1--;
+            if (xs<0)
+                xs=xs1=0;
 
             unsigned a = source(xs,ys);
             unsigned b = source(xs1,ys);
diff --git a/plugins/input/gdal/gdal_featureset.cpp b/plugins/input/gdal/gdal_featureset.cpp
index 768592bd0..05d30dae5 100644
--- a/plugins/input/gdal/gdal_featureset.cpp
+++ b/plugins/input/gdal/gdal_featureset.cpp
@@ -108,25 +108,30 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
     box2d<double> intersect = raster_extent.intersect(q.get_bbox());
     box2d<double> box = t.forward(intersect);
 
-    // TODO: error check this further...
-    float x_off_f = (intersect.minx()-raster_extent.minx()) / fabs(dx);
-    float y_off_f = (raster_extent.maxy()-intersect.maxy()) / fabs(dy);
-    
-    if (x_off_f < 0)
-    {
-	x_off_f = 0;
-    }
-    
-    if (y_off_f < 0)
-    {
-	y_off_f = 0;
-    }
-    
-    int x_off = static_cast<int>(x_off_f + 0.5);
-    int y_off = static_cast<int>(y_off_f + 0.5);
-   
-    int width = int(box.maxx() + 0.5) - int(box.minx() + 0.5);
-    int height = int(box.maxy() + 0.5) - int(box.miny() + 0.5);
+    //select minimum raster containing whole box
+    int x_off = static_cast<int>(floor(box.minx()));
+    int y_off = static_cast<int>(floor(box.miny()));
+    int end_x = static_cast<int>(ceil(box.maxx()));
+    int end_y = static_cast<int>(ceil(box.maxy()));
+    //clip to available data
+    if (x_off < 0)
+        x_off = 0;
+    if (y_off < 0)
+        y_off = 0;
+    if (end_x > raster_width)
+        end_x = raster_width;
+    if (end_y > raster_height)
+        end_y = raster_height;
+    int width = end_x - x_off;
+    int height = end_y - y_off;
+    // don't process almost invisible data
+    if (box.width() < 0.5)
+        width = 0;
+    if (box.height() < 0.5)
+        height = 0;
+    //calculate actual box2d of returned raster
+    box2d<double> feature_raster_extent(x_off, y_off, x_off+width, y_off+height); 
+    intersect = t.backward(feature_raster_extent);
     
 #ifdef MAPNIK_DEBUG         
     std::clog << "GDAL Plugin: Raster extent=" << raster_extent << "\n";
diff --git a/plugins/input/raster/raster_featureset.cpp b/plugins/input/raster/raster_featureset.cpp
index c03157225..1d8e9422d 100644
--- a/plugins/input/raster/raster_featureset.cpp
+++ b/plugins/input/raster/raster_featureset.cpp
@@ -73,21 +73,33 @@ feature_ptr raster_featureset<LookupPolicy>::next()
                CoordTransform t(image_width,image_height,extent_,0,0);
                box2d<double> intersect=bbox_.intersect(curIter_->envelope());
                box2d<double> ext=t.forward(intersect);
-               
-	       int start_x = int(ext.minx() + 0.5);
-	       int start_y = int(ext.miny() + 0.5);
-	       int end_x = int(ext.maxx() + 0.5);
-	       int end_y = int(ext.maxy() + 0.5);
-	       
-	       unsigned w = end_x - start_x;
-	       unsigned h = end_y - start_y;
-	       if ( w > 0 && h > 0)
-	       {
-		   image_data_32 image(w,h);
-		   reader->read(start_x,start_y,image);
-		   feature->set_raster(mapnik::raster_ptr(new raster(intersect,image)));
-	       }
-	    }
+               if ( ext.width()>0.5 && ext.height()>0.5 )
+               {
+                  //select minimum raster containing whole ext
+                  int x_off = static_cast<int>(floor(ext.minx()));
+                  int y_off = static_cast<int>(floor(ext.miny()));
+                  int end_x = static_cast<int>(ceil(ext.maxx()));
+                  int end_y = static_cast<int>(ceil(ext.maxy()));
+                  //clip to available data
+                  if (x_off < 0)
+                     x_off = 0;
+                  if (y_off < 0)
+                     y_off = 0;
+                  if (end_x > image_width)
+                     end_x = image_width;
+                  if (end_y > image_height)
+                     end_y = image_height;
+                  int width = end_x - x_off;
+                  int height = end_y - y_off;
+                  //calculate actual box2d of returned raster
+                  box2d<double> feature_raster_extent(x_off, y_off, x_off+width, y_off+height); 
+                  intersect = t.backward(feature_raster_extent);
+
+                  image_data_32 image(width,height);
+                  reader->read(x_off,y_off,image);
+                  feature->set_raster(mapnik::raster_ptr(new raster(intersect,image)));
+               }
+            }
          }
       }
       catch (...)
diff --git a/src/agg_renderer.cpp b/src/agg_renderer.cpp
index 8fe8e9175..7f25dd42b 100644
--- a/src/agg_renderer.cpp
+++ b/src/agg_renderer.cpp
@@ -738,12 +738,14 @@ void agg_renderer<T>::process(raster_symbolizer const& sym,
 	
 	box2d<double> ext=t_.forward(raster->ext_);
 	
-	int start_x = int(ext.minx() + 0.5);
-	int start_y = int(ext.miny() + 0.5);
-	int end_x = int(ext.maxx() + 0.5);
-	int end_y = int(ext.maxy() + 0.5);
-	int raster_width = end_x - start_x;
-	int raster_height = end_y - start_y;
+	int start_x = int(round(ext.minx()));
+	int start_y = int(round(ext.miny()));
+	int raster_width = int(round(ext.width()));
+	int raster_height = int(round(ext.height()));
+	int end_x = start_x + raster_width;
+	int end_y = start_y + raster_height;
+	double err_offs_x = (ext.minx()-start_x + ext.maxx()-end_x)/2;
+	double err_offs_y = (ext.miny()-start_y + ext.maxy()-end_y)/2;
 	
 	if ( raster_width > 0 && raster_height > 0)
 	{
@@ -752,9 +754,9 @@ void agg_renderer<T>::process(raster_symbolizer const& sym,
 	    if (sym.get_scaling() == "fast") {
 		scale_image<image_data_32>(target,raster->data_);
 	    } else if (sym.get_scaling() == "bilinear"){
-		scale_image_bilinear<image_data_32>(target,raster->data_);
+		scale_image_bilinear<image_data_32>(target,raster->data_, err_offs_x, err_offs_y);
 	    } else if (sym.get_scaling() == "bilinear8"){
-		scale_image_bilinear8<image_data_32>(target,raster->data_);
+		scale_image_bilinear8<image_data_32>(target,raster->data_, err_offs_x, err_offs_y);
 	    } else {
 		scale_image<image_data_32>(target,raster->data_);
 	    }
diff --git a/src/cairo_renderer.cpp b/src/cairo_renderer.cpp
index 7fb19eafa..d6aa65772 100644
--- a/src/cairo_renderer.cpp
+++ b/src/cairo_renderer.cpp
@@ -989,23 +989,33 @@ void cairo_renderer_base::process(raster_symbolizer const& sym,
     raster_ptr const& raster = feature.get_raster();
     if (raster)
     {
-	box2d<double> ext = t_.forward(raster->ext_);
-	image_data_32 target(int(ext.width() + 0.5), int(ext.height() + 0.5));
-	//TODO -- use cairo matrix transformation for scaling
-	if (sym.get_scaling() == "fast"){
-	    scale_image<image_data_32>(target, raster->data_);
-	} else if (sym.get_scaling() == "bilinear"){
-            scale_image_bilinear<image_data_32>(target,raster->data_);
-	} else if (sym.get_scaling() == "bilinear8"){
-            scale_image_bilinear8<image_data_32>(target,raster->data_);
-	} else {
-            scale_image<image_data_32>(target,raster->data_);
-	}
+        box2d<double> ext = t_.forward(raster->ext_);
+        int start_x = int(round(ext.minx()));
+        int start_y = int(round(ext.miny()));
+        int raster_width = int(round(ext.width()));
+        int raster_height = int(round(ext.height()));
+        int end_x = start_x + raster_width;
+        int end_y = start_y + raster_height;
+        double err_offs_x = (ext.minx()-start_x + ext.maxx()-end_x)/2;
+        double err_offs_y = (ext.miny()-start_y + ext.maxy()-end_y)/2;
 
-	cairo_context context(context_);
-
-	//TODO -- support for advanced image merging
-	context.add_image(int(ext.minx()+0.5), int(ext.miny()+0.5), target, sym.get_opacity());
+        if (raster_width > 0 && raster_height > 0)
+        {
+            image_data_32 target(raster_width, raster_height);
+            //TODO -- use cairo matrix transformation for scaling
+            if (sym.get_scaling() == "fast"){
+                scale_image<image_data_32>(target, raster->data_);
+            } else if (sym.get_scaling() == "bilinear"){
+                scale_image_bilinear<image_data_32>(target,raster->data_, err_offs_x, err_offs_y);
+            } else if (sym.get_scaling() == "bilinear8"){
+                scale_image_bilinear8<image_data_32>(target,raster->data_, err_offs_x, err_offs_y);
+            } else {
+                scale_image<image_data_32>(target,raster->data_);
+            }
+            cairo_context context(context_);
+            //TODO -- support for advanced image merging
+            context.add_image(start_x, start_y, target, sym.get_opacity());
+        }
     }
 }