A fix for two distinct issues associated with gdal featuresets, the first is overviews were not properly being utilized based on the resolution of the final image requested. The second is that allocation of far too much memory could be possible in GDAL to attempt to do resampling internally in mapnik. This now has a hard cap so that we allocate less memory in these situations but are still able to resample internally in mapnik.

include/mapnik/image.hpp

@@ -74,7 +74,7 @@ public:
     static constexpr image_dtype dtype = T::id;
     static constexpr std::size_t pixel_size = sizeof(pixel_type);
 private:
-    detail::image_dimensions<65535> dimensions_;
+    detail::image_dimensions<4294836225> dimensions_;
     detail::buffer buffer_;
     double offset_;
     double scaling_;

include/mapnik/image_impl.hpp

@@ -39,8 +39,10 @@ image_dimensions<max_size>::image_dimensions(int width, int height)
     : width_(width),
       height_(height)
 {
-    if (width < 0 || static_cast<std::size_t>(width) > max_size) throw std::runtime_error("Invalid width for image dimensions requested");
+    int64_t area = (int64_t)width * (int64_t)height;
-    if (height < 0 || static_cast<std::size_t>(height) > max_size) throw std::runtime_error("Invalid height for image dimensions requested");
+    if (width < 0) throw std::runtime_error("Invalid width for image dimensions requested");
+    if (height < 0) throw std::runtime_error("Invalid height for image dimensions requested");
+    if (area > max_size) throw std::runtime_error("Image area too large based on image dimensions");
 }
 
 template <std::size_t max_size>

plugins/input/gdal/gdal_featureset.cpp

@@ -155,6 +155,8 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
     //size of resized output pixel in source image domain
     double margin_x = 1.0 / (std::fabs(dx_) * std::get<0>(q.resolution()));
     double margin_y = 1.0 / (std::fabs(dy_) * std::get<1>(q.resolution()));
+    MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: margin_x=" << margin_x;
+    MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: margin_y=" << margin_y;
     if (margin_x < 1)
     {
         margin_x = 1.0;
@@ -169,6 +171,10 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
     int y_off = rint(box.miny() - margin_y);
     int end_x = rint(box.maxx() + margin_x);
     int end_y = rint(box.maxy() + margin_y);
+    MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: x_off=" << x_off;
+    MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: y_off=" << y_off;
+    MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: end_x=" << end_x;
+    MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: end_y=" << end_y;
 
     //clip to available data
     if (x_off < 0)
@@ -187,21 +193,139 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
     {
         end_y = raster_height_;
     }
+
+    // width and height of the portion of the source image we are requesting
     int width = end_x - x_off;
     int height = end_y - y_off;
 
+    // In many cases we want GDAL to simply return the exact image so we
+    // can handle resampling internally in mapnik. In other cases such as 
+    // when overviews exist or when the image allocated might be too large
+    // we want to utilize some resampling in GDAL instead.
+    int im_height = height;
+    int im_width = width;
+    double im_offset_x = x_off;
+    double im_offset_y = y_off;
+    int current_width = (int)raster_width_;
+    int current_height = (int)raster_height_;
+
+    // loop through overviews -- snap up in resolution to closest overview if necessary
+    // we find an image size that most resembles the resolution of our output image.
+    double width_res = std::get<0>(q.resolution());
+    double height_res = std::get<1>(q.resolution());
+    int res_adjusted_raster_width = static_cast<int>(std::floor(((double)raster_width_ * width_res) + .5));
+    int res_adjusted_raster_height = static_cast<int>(std::floor(((double)raster_height_ * height_res) + .5));
+    if (band_ > 0 && band_ < nbands_)
+    {
+        GDALRasterBand * band = dataset_.GetRasterBand(band_);
+        int band_overviews = band->GetOverviewCount();
+        if (band_overviews > 0)
+        {
+            for (int b = 0; b < band_overviews; b++)
+            {
+                GDALRasterBand * overview = band->GetOverview(b);
+                int overview_width = overview->GetXSize();
+                int overview_height = overview->GetYSize();
+                if ((overview_width < current_width || overview_height < current_height) &&
+                    res_adjusted_raster_width <= overview_width &&
+                    res_adjusted_raster_height <= overview_height)
+                {
+                    current_width = overview_width;
+                    current_height = overview_height;
+                }
+            }
+        }
+    }
+    else
+    {
+        for (int i = 0; i < nbands_; ++i)
+        {
+            GDALRasterBand * band = dataset_.GetRasterBand(i + 1);
+            int band_overviews = band->GetOverviewCount();
+            if (band_overviews > 0)
+            {
+                for (int b = 0; b < band_overviews; b++)
+                {
+                    GDALRasterBand * overview = band->GetOverview(b);
+                    int overview_width = overview->GetXSize();
+                    int overview_height = overview->GetYSize();
+                    if ((overview_width < current_width || overview_height < current_height) &&
+                        res_adjusted_raster_width <= overview_width &&
+                        res_adjusted_raster_height <= overview_height)
+                    {
+                        current_width = overview_width;
+                        current_height = overview_height;
+                    }
+                }
+            }
+        }
+    }
+    if (current_width != (int)raster_width_ || current_height != (int)raster_height_)
+    {
+        if (current_width != (int)raster_width_)
+        {
+            double ratio = (double)current_width / (double)raster_width_;
+            int adjusted_width = static_cast<int>(std::floor((ratio * im_width) + 0.5));
+            double adjusted_ratio = (double)adjusted_width / (double)im_width;
+            im_offset_x = adjusted_ratio * im_offset_x;
+            im_width = adjusted_width;
+        }
+        if (current_height != (int)raster_height_)
+        {
+            double ratio = (double)current_height / (double)raster_height_;
+            int adjusted_height = static_cast<int>(std::floor((ratio * im_height) + 0.5));
+            double adjusted_ratio = (double)adjusted_height / (double)im_height;
+            im_offset_y = adjusted_ratio * im_offset_y;
+            im_height = adjusted_height;
+        }
+    }
+    
+    // Maximum memory limitation for image will be simply based on the maximum
+    // area we allow for an image. The true memory footprint therefore will vary based
+    // on the type of imagery that exists. 
+    // max_im_area based on 50 mb limit for RGBA
+    constexpr int64_t max_im_area = (50*1024*1024) / 4;
+    int64_t im_area = (int64_t)im_width * (int64_t)im_height;
+    if (im_area > max_im_area)
+    {
+        int adjusted_width = static_cast<int>(std::round(std::sqrt(max_im_area * ((double)im_width / (double)im_height))));
+        int adjusted_height = static_cast<int>(std::round(std::sqrt(max_im_area * ((double)im_height / (double)im_width))));
+        if (adjusted_width < 1)
+        {
+            adjusted_width = 1;
+        }
+        if (adjusted_height < 1)
+        {
+            adjusted_height = 1;
+        }
+        double ratio_x = (double)adjusted_width / (double)im_width;
+        double ratio_y = (double)adjusted_height / (double)im_height;
+        im_offset_x = ratio_x * im_offset_x;
+        im_offset_y = ratio_y * im_offset_y;
+        im_width = adjusted_width;
+        im_height = adjusted_height;
+        current_width = static_cast<int>(std::floor((ratio_x * current_width) + 0.5));
+        current_height = static_cast<int>(std::floor((ratio_y * current_height) + 0.5));
+    }
+    
     //calculate actual box2d of returned raster
-    box2d<double> feature_raster_extent(x_off, y_off, x_off + width, y_off + height);
+    view_transform t2(current_width, current_height, raster_extent_, 0, 0);
-    feature_raster_extent = t.backward(feature_raster_extent);
+    box2d<double> feature_raster_extent(im_offset_x, im_offset_y, im_offset_x + im_width, im_offset_y + im_height);
+    MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Feature Raster extent=" << feature_raster_extent;
+    feature_raster_extent = t2.backward(feature_raster_extent);
 
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Raster extent=" << raster_extent_;
+    MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Feature Raster extent=" << feature_raster_extent;
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: View extent=" << intersect;
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Query resolution=" << std::get<0>(q.resolution()) << "," << std::get<1>(q.resolution());
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: StartX=" << x_off << " StartY=" << y_off << " Width=" << width << " Height=" << height;
+    MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: IM StartX=" << im_offset_x << " StartY=" << im_offset_y << " Width=" << im_width << " Height=" << im_height;
 
     if (width > 0 && height > 0)
     {
-        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Image Size=(" << width << "," << height << ")";
+
+        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Requested Image Size=(" << width << "," << height << ")";
+        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Image Size=(" << im_width << "," << im_height << ")";
         MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Reading band=" << band_;
         if (band_ > 0) // we are querying a single band
         {
@@ -217,7 +341,7 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
             {
             case GDT_Byte:
             {
-                mapnik::image_gray8 image(width, height);
+                mapnik::image_gray8 image(im_width, im_height);
                 image.set(std::numeric_limits<std::uint8_t>::max());
                 raster_nodata = band->GetNoDataValue(&raster_has_nodata);
                 raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
@@ -237,7 +361,7 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
             case GDT_Float64:
             case GDT_Float32:
             {
-                mapnik::image_gray32f image(width, height);
+                mapnik::image_gray32f image(im_width, im_height);
                 image.set(std::numeric_limits<float>::max());
                 raster_nodata = band->GetNoDataValue(&raster_has_nodata);
                 raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
@@ -256,7 +380,7 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
             }
             case GDT_UInt16:
             {
-                mapnik::image_gray16 image(width, height);
+                mapnik::image_gray16 image(im_width, im_height);
                 image.set(std::numeric_limits<std::uint16_t>::max());
                 raster_nodata = band->GetNoDataValue(&raster_has_nodata);
                 raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
@@ -276,7 +400,7 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
             default:
             case GDT_Int16:
             {
-                mapnik::image_gray16s image(width, height);
+                mapnik::image_gray16s image(im_width, im_height);
                 image.set(std::numeric_limits<std::int16_t>::max());
                 raster_nodata = band->GetNoDataValue(&raster_has_nodata);
                 raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
@@ -297,7 +421,7 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
         }
         else // working with all bands
         {
-            mapnik::image_rgba8 image(width, height);
+            mapnik::image_rgba8 image(im_width, im_height);
             image.set(std::numeric_limits<std::uint32_t>::max());
             for (int i = 0; i < nbands_; ++i)
             {

src/image.cpp

@@ -84,7 +84,7 @@ void buffer::swap(buffer & rhs)
     std::swap(owns_, rhs.owns_);
 }
 
-template struct MAPNIK_DECL image_dimensions<65535>;
+template struct MAPNIK_DECL image_dimensions<4294836225>;
 
 } // end ns detail
 

test/data

@@ -1 +1 @@
-Subproject commit da98ee10eb7ac8cd6d3b1441c3470dc5d3748c6e
+Subproject commit bafeaab6922051a5114648e5ccffe4104e23081d

test/data-visual

@@ -1 +1 @@
-Subproject commit 9e0887293e7f6195e1eecc6fed13ecb64e300259
+Subproject commit b25e93ad4cd4e29b07384e4a6949dff506c3d939