On Sat, 2008-10-25 at 22:56 +0100, Jon Burgess wrote:
> Traceback (most recent call last):
>   File "./nik2img.py", line 571, in <module>
>     render(mapnik_map,o,FORMAT)
>   File "./nik2img.py", line 188, in render
>     mapnik.render_to_file(*args)
> MemoryError
> 
> I'm not certain at this point whether this is a bug in the test setup
> or
> a fundamental flaw in the Mapnik GDAL handling. I guess the approach
> of
> turning the source file into a bitmap works fine for normal sized
> images
> but falls down when faced with such an enormous raw bitmap.
> 
> Mapnik would be much better off it could extract lower resolution
> image
> from the multi-resolution file. I don't know if GDAL exports this
> aspect
> of the file through the API. If not, it might need a special MrSID
> input
> plugin in Mapnik.
> 

GDAL exposes the lower resolutions via overviews. I put together a quick
test and it appears to work. The patch is attached. It is definitely
_not_ ready to be merged. Currently just picks the smallest overview
greater than 1000 pixels wide. 

In reality it needs to be selecting the best overview for the current
query and I've no idea how this will work for other GDAL image formats.

        Jon


Index: plugins/input/gdal/gdal_featureset.cpp
===================================================================
--- plugins/input/gdal/gdal_featureset.cpp	(revision 750)
+++ plugins/input/gdal/gdal_featureset.cpp	(working copy)
@@ -30,10 +30,11 @@
 using mapnik::feature_ptr;
 using mapnik::CoordTransform;
 
-gdal_featureset::gdal_featureset(GDALDataset & dataset, query const& q)
+gdal_featureset::gdal_featureset(GDALDataset & dataset, query const& q, int overview)
    : dataset_(dataset),
      query_extent_(q.get_bbox()),
-     first_(true) {}
+     first_(true),
+     overview_(overview) {}
 
 gdal_featureset::~gdal_featureset() {}
 
@@ -48,14 +49,20 @@
       GDALRasterBand * blue = 0;
       GDALRasterBand * alpha = 0;
       GDALRasterBand * grey = 0; 
+      unsigned raster_xsize = 0;
+      unsigned raster_ysize = 0;
+
       for (int i = 0; i < dataset_.GetRasterCount() ;++i)
       {
          GDALRasterBand * band = dataset_.GetRasterBand(i+1);
-         //if (band->GetOverviewCount() > 0)
-         //{
-         //  band = band->GetOverview(0);
-         //}
-         
+         if (overview_ >= 0)
+	   band = band->GetOverview(overview_);
+
+	 if (raster_xsize == 0)
+	   raster_xsize =  band->GetXSize();
+	 if (raster_ysize == 0)
+	   raster_ysize =  band->GetYSize();
+
 #ifdef MAPNIK_DEBUG         
          int bsx,bsy;
          band->GetBlockSize(&bsx,&bsy);
@@ -112,8 +119,8 @@
          }
       }
       
-      unsigned raster_xsize = dataset_.GetRasterXSize();
-      unsigned raster_ysize = dataset_.GetRasterYSize();
+      //unsigned raster_xsize = dataset_.GetRasterXSize();
+      //unsigned raster_ysize = dataset_.GetRasterYSize();
       double tr[6];
       dataset_.GetGeoTransform(tr);
       double x0 = tr[0];
@@ -125,6 +132,14 @@
       Envelope<double> intersection = raster_extent.intersect(query_extent_);
       Envelope<double> box = t.forward(intersection);
       
+      std::cerr << "GDAL gdal_featureset::next(): Raster Image(" << raster_xsize << "," << raster_ysize << ")\n";
+
+      std::cerr << "GDAL gdal_featureset::next(): Raster " << raster_extent << "\n";
+      std::cerr << "GDAL gdal_featureset::next(): Query " << query_extent_ << "\n";
+
+      std::cerr << "GDAL gdal_featureset::next(): Raster Extent(" << raster_extent.width() << "," << raster_extent.height() << ")\n";
+      std::cerr << "GDAL gdal_featureset::next(): Query Extent(" << query_extent_.width() << "," << query_extent_.height() << ")\n";
+
       int start_x = int(box.minx());
       int start_y = int(box.miny());
       int width = int(box.width());
@@ -132,6 +147,7 @@
       
       if (width > 0 && height > 0)
       {
+	std::cerr << "GDAL gdal_featureset::next(): Creating image(" << width << "," << height << ")\n";
          mapnik::ImageData32 image(width,height);
          image.set(0xffffffff);
          
Index: plugins/input/gdal/gdal_datasource.cpp
===================================================================
--- plugins/input/gdal/gdal_datasource.cpp	(revision 750)
+++ plugins/input/gdal/gdal_datasource.cpp	(working copy)
@@ -47,12 +47,35 @@
 
    dataset_ = boost::shared_ptr<GDALDataset>(reinterpret_cast<GDALDataset*>(GDALOpen((*file).c_str(),GA_ReadOnly)));
    if (!dataset_) throw datasource_exception("failed to create GDALDataset");
+
+   // FIXME: Down res to some sensible size
+   overview_ = -1;
+   GDALRasterBand * band = dataset_->GetRasterBand(1);
+   if (band->GetOverviewCount() > 0) {
+     int j;
+     std::cerr << "Overview count = " << band->GetOverviewCount() << "\n";
+     for (j = 0; j < band->GetOverviewCount(); j++) {
+       GDALRasterBand * oview = band->GetOverview(j); 
+       std::cerr << "Size: " << oview->GetXSize() << "," << oview->GetXSize() << "\n";
+       if (oview->GetXSize() < 1000)
+	 break;
+     }
+     overview_ = j-1;
+     band = band->GetOverview(overview_);
+     std::cerr << "Picked overview " << overview_ << " (" << band->GetXSize() << "," << band->GetXSize() << ")\n";
+   }
+
+   xsize_ = band->GetXSize();
+   ysize_ = band->GetYSize();
+
    double tr[6];
    dataset_->GetGeoTransform(tr);
    double x0 = tr[0];
    double y0 = tr[3];
-   double x1 = tr[0] + dataset_->GetRasterXSize()*tr[1] + dataset_->GetRasterYSize()*tr[2];
-   double y1 = tr[3] + dataset_->GetRasterXSize()*tr[4] + dataset_->GetRasterYSize()*tr[5];
+   //double x1 = tr[0] + dataset_->GetRasterXSize()*tr[1] + dataset_->GetRasterYSize()*tr[2];
+   //double y1 = tr[3] + dataset_->GetRasterXSize()*tr[4] + dataset_->GetRasterYSize()*tr[5];
+   double x1 = tr[0] + xsize_*tr[1] + ysize_*tr[2];
+   double y1 = tr[3] + xsize_*tr[4] + ysize_*tr[5];
    extent_.init(x0,y0,x1,y1);
 }
 
@@ -82,7 +105,7 @@
 {
    if (dataset_)
    {
-      return featureset_ptr(new gdal_featureset(*dataset_, q));
+     return featureset_ptr(new gdal_featureset(*dataset_, q, overview_));
    }
    return featureset_ptr();
 }
Index: plugins/input/gdal/gdal_featureset.hpp
===================================================================
--- plugins/input/gdal/gdal_featureset.hpp	(revision 750)
+++ plugins/input/gdal/gdal_featureset.hpp	(working copy)
@@ -32,13 +32,14 @@
 {
    public:
       
-      gdal_featureset(GDALDataset & dataset,mapnik::query const& q);
+      gdal_featureset(GDALDataset & dataset,mapnik::query const& q, int overview);
       virtual ~gdal_featureset();
       mapnik::feature_ptr next();
    private:
       GDALDataset & dataset_;
       mapnik::Envelope<double> query_extent_;
       bool first_;
+      int overview_;
 };
 
 #endif // GDAL_FEATURESET_HPP
Index: plugins/input/gdal/gdal_datasource.hpp
===================================================================
--- plugins/input/gdal/gdal_datasource.hpp	(revision 750)
+++ plugins/input/gdal/gdal_datasource.hpp	(working copy)
@@ -43,6 +43,8 @@
       mapnik::Envelope<double> extent_;
       boost::shared_ptr<GDALDataset> dataset_;
       mapnik::layer_descriptor desc_;
+      int overview_;
+      int xsize_, ysize_;
 };
 
 
_______________________________________________
Mapnik-users mailing list
[email protected]
https://lists.berlios.de/mailman/listinfo/mapnik-users

Reply via email to