diff --git a/plugins/input/gdal/gdal_featureset.cpp b/plugins/input/gdal/gdal_featureset.cpp
index 2f9dbbbcd..55bc61a9f 100644
--- a/plugins/input/gdal/gdal_featureset.cpp
+++ b/plugins/input/gdal/gdal_featureset.cpp
@@ -33,6 +33,7 @@ using mapnik::CoordTransform;
 using mapnik::point_impl;
 using mapnik::geometry2d;
 
+#define MAPNIK_DEBUG
 
 gdal_featureset::gdal_featureset(GDALDataset & dataset, int band, gdal_query q)
    : dataset_(dataset),
@@ -43,7 +44,7 @@ gdal_featureset::gdal_featureset(GDALDataset & dataset, int band, gdal_query q)
 gdal_featureset::~gdal_featureset()
 {
 #ifdef MAPNIK_DEBUG
-      std::clog << "closing dataset = " << &dataset_ << "\n";
+      std::clog << "GDAL Plugin: closing dataset = " << &dataset_ << "\n";
 #endif
    GDALClose(&dataset_);
 }
@@ -54,7 +55,7 @@ feature_ptr gdal_featureset::next()
    {
       first_ = false;
 #ifdef MAPNIK_DEBUG
-      std::clog << "featureset, dataset = " << &dataset_ << "\n";
+      std::clog << "GDAL Plugin: featureset, dataset = " << &dataset_ << "\n";
 #endif
 
       query *q = boost::get<query>(&gquery_);
@@ -71,7 +72,6 @@ feature_ptr gdal_featureset::next()
    return feature_ptr();
 }
 
-       
 feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
 {
    feature_ptr feature(new Feature(1));
@@ -82,98 +82,85 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
    GDALRasterBand * alpha = 0;
    GDALRasterBand * grey = 0; 
 
-   unsigned raster_xsize = dataset_.GetRasterXSize();
-   unsigned raster_ysize = dataset_.GetRasterYSize();
    int nbands = dataset_.GetRasterCount();
-
+   
+   unsigned raster_width = dataset_.GetRasterXSize();
+   unsigned raster_height = dataset_.GetRasterYSize();
+   
+   // TODO - pull from class attributes...
    double tr[6];
    dataset_.GetGeoTransform(tr);
    double dx = tr[1];
    double dy = tr[5];
-   double x0 = tr[0];
-   double y0 = tr[3];
-   double x1 = tr[0] + raster_xsize * dx + raster_ysize * tr[2];
-   double y1 = tr[3] + raster_xsize * tr[4] + raster_ysize * dy;
-   Envelope<double> raster_extent(x0,y0,x1,y1);
-   CoordTransform t (raster_xsize,raster_ysize,raster_extent,0,0);
-   Envelope<double> intersection = raster_extent.intersect(q.get_bbox());
-   Envelope<double> box = t.forward(intersection);
+   double x0 = tr[0]; // minx
+   double y0 = tr[3]; // miny
+   double x1 = tr[0] + raster_width * dx + raster_height * tr[2]; // maxx
+   double y1 = tr[3] + raster_width * tr[4] + raster_height * dy; // maxy
+   Envelope<double> raster_extent(x0,y0,x1,y1); 
+   
+   CoordTransform t(raster_width,raster_height,raster_extent,0,0);
+   Envelope<double> intersect = raster_extent.intersect(q.get_bbox());
+   Envelope<double> box = t.forward(intersect);
 
-   int start_x, start_y;
-   if (dx>0)
-       start_x = int(std::floor(intersection.minx()-raster_extent.minx())/dx);
-   else
-       // west-left raster
-       start_x = -int(std::floor(raster_extent.maxx()-intersection.maxx())/dx);
-   if(dy>0)
-       start_y = int(std::floor(intersection.miny()-raster_extent.miny())/dy);
-   else
-       // north-up raster
-       start_y = -int(std::floor(raster_extent.maxy()-intersection.maxy())/dy);
-
-   int width = int(box.width()+0.5);
-   int height = int(box.height()+0.5);
+   // 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);
+   int y_off = static_cast<int>(y_off_f);
 
+   int width = int(box.width() + 0.5);
+   int height = int(box.height() + 0.5);
+    
 #ifdef MAPNIK_DEBUG         
-   std::cout << "raster_extent=" << raster_extent << "\n";
-   std::cout << "intersection=" << intersection << "\n";
-   std::cout << boost::format("RasterWidth=%d RasterHeight=%d \n") % raster_xsize % raster_ysize;
-   std::cout << boost::format("StartX=%d StartY=%d Width=%d Height=%d \n") % start_x % start_y % width % height;
+   std::cout << "GDAL Plugin: Raster extent=" << raster_extent << "\n";
+   std::cout << "GDAL Plugin: View extent=" << intersect << "\n";
+   std::cout << "GDAL Plugin: Query resolution=" << q.resolution() << "\n";
+   std::cout << boost::format("GDAL Plugin: StartX=%d StartY=%d Width=%d Height=%d \n") % x_off % y_off % width % height;
 #endif
 
    if (width > 0 && height > 0)
    {
-      int im_width = q.resolution()*intersection.width();
-      int im_height = q.resolution()*intersection.height();
+      int im_width = int(q.resolution() * intersect.width() + 0.5);
+      int im_height = int(q.resolution() * intersect.height() + 0.5);
 
       mapnik::ImageData32 image(im_width, im_height);
       image.set(0xffffffff);
          
-#ifdef MAPNIK_DEBUG         
-      std::cout << "band=" << band_ << "\n";
-      std::cout << "im_size=(" << im_width << "," << im_height << ")\n";
-      std::cout << "box_size=(" << width << "," << height << ")\n";
+#ifdef MAPNIK_DEBUG
+      std::cout << "GDAL Plugin: Image Size=(" << im_width << "," << im_height << ")\n";
+      std::cout << "GDAL Plugin: Reading band " << band_ << "\n";
 #endif
 
-      if(band_>0)      // we are querying a single band
+      if (band_>0) // we are querying a single band
       {
          float *imageData = (float*)image.getBytes();
          GDALRasterBand * band = dataset_.GetRasterBand(band_);
-         band->RasterIO(GF_Read, start_x, start_y, width, height,
+         band->RasterIO(GF_Read, x_off, y_off, width, height,
                         imageData, image.width(), image.height(),
                         GDT_Float32, 0, 0);
 
-         feature->set_raster(mapnik::raster_ptr(new mapnik::raster(intersection,image)));
+         feature->set_raster(mapnik::raster_ptr(new mapnik::raster(intersect,image)));
       }
       
-      else             // working with all bands
+      else // working with all bands
       {
          for (int i = 0; i < nbands; ++i)
          {
             GDALRasterBand * band = dataset_.GetRasterBand(i+1);
          
 #ifdef MAPNIK_DEBUG
-            int band_overviews = band->GetOverviewCount();
-            // Note: overviews are picked automatically during GdalDataset::RasterIO below
-            // when nXSize > nPixelSpace and nYSize > nLineSpace 
-            
-            if (band_overviews > 0)
-            {
-                for (int b = 0; b < band_overviews; b++)
-                {
-                    GDALRasterBand * overview = band->GetOverview (b);
-                    std::cout << boost::format("Overview=%d Width=%d Height=%d \n")
-                                           % b % overview->GetXSize() % overview->GetYSize();
-                }
-            }
-
-            int bsx,bsy;
-            double scale;
-            band->GetBlockSize(&bsx,&bsy);
-            scale = band->GetScale();
-            std::cout << boost::format("Block=%dx%d Scale=%f Type=%s Color=%s \n") % bsx % bsy % scale
-               % GDALGetDataTypeName(band->GetRasterDataType())
-               % GDALGetColorInterpretationName(band->GetColorInterpretation());   
+            get_overview_meta(band);  
 #endif
 
             GDALColorInterp color_interp = band->GetColorInterpretation();
@@ -181,68 +168,108 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
             {
                case GCI_RedBand:
                   red = band;
+#ifdef MAPNIK_DEBUG
+                  std::cout << "GDAL Plugin: Found red band" << "\n";
+#endif
                   break;
                case GCI_GreenBand:
                   green = band;
+#ifdef MAPNIK_DEBUG
+                  std::cout << "GDAL Plugin: Found green band" << "\n";
+#endif
                   break;
                case GCI_BlueBand:
                   blue = band;
+#ifdef MAPNIK_DEBUG
+                  std::cout << "GDAL Plugin: Found blue band" << "\n";
+#endif
                   break;
                case GCI_AlphaBand:
                   alpha = band;
+#ifdef MAPNIK_DEBUG
+                  std::cout << "GDAL Plugin: Found alpha band" << "\n";
+#endif
                   break;
                case GCI_GrayIndex:
                   grey = band;
+#ifdef MAPNIK_DEBUG
+                  std::cout << "GDAL Plugin: Found gray band" << "\n";
+#endif
                   break;
                case GCI_PaletteIndex:
                {
                   grey = band;
+#ifdef MAPNIK_DEBUG
+                  std::cout << "GDAL Plugin: Found gray band, and colortable..." << "\n";
+#endif
                   GDALColorTable *color_table = band->GetColorTable();
             
                   if ( color_table)
                   {
                      int count = color_table->GetColorEntryCount();
 #ifdef MAPNIK_DEBUG
-                     std::cout << "Color Table count = " << count << "\n";
+                     std::cout << "GDAL Plugin: Color Table count = " << count << "\n";
 #endif 
                      for ( int i = 0; i < count; i++ )
                      {
                         const GDALColorEntry *ce = color_table->GetColorEntry ( i );
                         if (!ce ) continue;
 #ifdef MAPNIK_DEBUG
-                        std::cout << "color entry RGB (" << ce->c1 <<"," <<ce->c2 << "," << ce->c3 << ")\n"; 
+                        std::cout << "GDAL Plugin: Color entry RGB (" << ce->c1 << "," <<ce->c2 << "," << ce->c3 << ")\n"; 
 #endif
                      }
                   }
                   break;
                }
                case GCI_Undefined:
+#ifdef MAPNIK_DEBUG
+                  std::cout << "GDAL Plugin: Found undefined band (assumming gray band)" << "\n";
+#endif
                   grey = band;
                   break;
                default:
+#ifdef MAPNIK_DEBUG
+                  std::cout << "GDAL Plugin: band type unknown!" << "\n";
+#endif
                   break;
             }
          }
 
-         if (red && green && blue)
+          if (red && green && blue)
          {
-            red->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 0, image.width(),image.height(), GDT_Byte,4,4*image.width());
-            green->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 1, image.width(),image.height(), GDT_Byte,4,4*image.width());
-            blue->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 2, image.width(),image.height(), GDT_Byte,4,4*image.width());
+#ifdef MAPNIK_DEBUG
+            std::cout << "GDAL Plugin: processing rgb bands..." << "\n";
+#endif
+            red->RasterIO(GF_Read,x_off,y_off,width,height,image.getBytes() + 0,
+                image.width(),image.height(),GDT_Byte,4,4*image.width());
+            green->RasterIO(GF_Read,x_off,y_off,width,height,image.getBytes() + 1,
+                image.width(),image.height(),GDT_Byte,4,4*image.width());
+            blue->RasterIO(GF_Read,x_off,y_off,width,height,image.getBytes() + 2,
+                image.width(),image.height(),GDT_Byte,4,4*image.width());
          }
          else if (grey)
          {
-            grey->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 0, image.width(),image.height(), GDT_Byte,4,4*image.width());
-            grey->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 1, image.width(),image.height(), GDT_Byte,4,4*image.width());
-            grey->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 2, image.width(),image.height(), GDT_Byte,4,4*image.width());
+#ifdef MAPNIK_DEBUG
+            std::cout << "GDAL Plugin: processing gray band..." << "\n";
+#endif
+            grey->RasterIO(GF_Read,x_off,y_off,width,height,image.getBytes() + 0,
+                image.width(),image.height(),GDT_Byte, 4, 4 * image.width());
+            grey->RasterIO(GF_Read,x_off,y_off,width,height,image.getBytes() + 1,
+                image.width(),image.height(),GDT_Byte, 4, 4 * image.width());
+            grey->RasterIO(GF_Read,x_off,y_off,width,height,image.getBytes() + 2,
+                image.width(),image.height(),GDT_Byte, 4, 4 * image.width());
          }
 
          if (alpha)
          {
-            alpha->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 3, image.width(),image.height(), GDT_Byte,4,4*image.width());
+#ifdef MAPNIK_DEBUG
+            std::cout << "GDAL Plugin: processing alpha band..." << "\n";
+#endif
+            alpha->RasterIO(GF_Read,x_off,y_off,width,height,image.getBytes() + 3,
+                image.width(),image.height(),GDT_Byte,4,4*image.width());
          }
 
-         feature->set_raster(mapnik::raster_ptr(new mapnik::raster(intersection,image)));
+         feature->set_raster(mapnik::raster_ptr(new mapnik::raster(intersect,image)));
       }
       return feature;
    }
@@ -250,7 +277,6 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
 }
 
 
-
 feature_ptr gdal_featureset::get_feature_at_point(mapnik::coord2d const& pt)
 {
    if (band_>0) {
@@ -272,8 +298,8 @@ feature_ptr gdal_featureset::get_feature_at_point(mapnik::coord2d const& pt)
 
       if(0<=x && x<raster_xsize && 0<=y && y<raster_ysize) {
 #ifdef MAPNIK_DEBUG         
-         std::cout << boost::format("pt.x=%f pt.y=%f\n") % pt.x % pt.y;
-         std::cout << boost::format("x=%f y=%f\n") % x % y;
+         std::cout << boost::format("GDAL Plugin: pt.x=%f pt.y=%f\n") % pt.x % pt.y;
+         std::cout << boost::format("GDAL Plugin: x=%f y=%f\n") % x % y;
 #endif
          GDALRasterBand * band = dataset_.GetRasterBand(band_);
          int hasNoData;
@@ -293,3 +319,30 @@ feature_ptr gdal_featureset::get_feature_at_point(mapnik::coord2d const& pt)
    }
    return feature_ptr();
 }
+
+void gdal_featureset::get_overview_meta(GDALRasterBand * band)
+{
+    int band_overviews = band->GetOverviewCount();    
+    if (band_overviews > 0)
+    {
+        std::cout << "GDAL Plugin: "<< band_overviews << " overviews found!" << "\n";
+        for (int b = 0; b < band_overviews; b++)
+        {
+            GDALRasterBand * overview = band->GetOverview (b);
+            std::cout << boost::format("GDAL Plugin: Overview=%d Width=%d Height=%d \n")
+                                   % b % overview->GetXSize() % overview->GetYSize();
+        }
+    }
+    else
+    {
+        std::cout << "GDAL Plugin: No overviews found!" << "\n";
+    }
+    
+    int bsx,bsy;
+    double scale;
+    band->GetBlockSize(&bsx,&bsy);
+    scale = band->GetScale();
+    std::cout << boost::format("GDAL Plugin: Block=%dx%d Scale=%f Type=%s Color=%s \n") % bsx % bsy % scale
+       % GDALGetDataTypeName(band->GetRasterDataType())
+       % GDALGetColorInterpretationName(band->GetColorInterpretation()); 
+}
diff --git a/plugins/input/gdal/gdal_featureset.hpp b/plugins/input/gdal/gdal_featureset.hpp
index 0e4d4abc6..ce0bfba48 100644
--- a/plugins/input/gdal/gdal_featureset.hpp
+++ b/plugins/input/gdal/gdal_featureset.hpp
@@ -28,6 +28,7 @@
 #include <boost/variant.hpp>
 
 class GDALDataset;
+class GDALRasterBand;
 
 typedef boost::variant<mapnik::query,mapnik::coord2d> gdal_query;
 
@@ -41,6 +42,7 @@ class gdal_featureset : public mapnik::Featureset
    private:
       mapnik::feature_ptr get_feature(mapnik::query const& q);
       mapnik::feature_ptr get_feature_at_point(mapnik::coord2d const& p);
+      void get_overview_meta(GDALRasterBand * band);
       GDALDataset & dataset_;
       int band_;
       gdal_query gquery_;