mapnik/plugins/input/gdal/gdal_featureset.cpp

/*****************************************************************************
 * 
 * This file is part of Mapnik (c++ mapping toolkit)
 *
 * Copyright (C) 2007 Artem Pavlenko
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 *****************************************************************************/
//$Id$

#include "gdal_featureset.hpp"
#include <gdal_priv.h>

// boost
#include <boost/format.hpp>

using mapnik::query;
using mapnik::coord2d;
using mapnik::Envelope;
using mapnik::Feature;
using mapnik::feature_ptr;
using mapnik::CoordTransform;
using mapnik::point_impl;
using mapnik::geometry2d;


gdal_featureset::gdal_featureset(GDALDataset & dataset, int band, gdal_query q)
   : dataset_(dataset),
     band_(band),
     gquery_(q),
     first_(true) {}

gdal_featureset::~gdal_featureset()
{
#ifdef MAPNIK_DEBUG
   std::clog << "GDAL Plugin: closing dataset = " << &dataset_ << "\n";
#endif
   GDALClose(&dataset_);
}

feature_ptr gdal_featureset::next()
{
   if (first_)
   {
      first_ = false;
#ifdef MAPNIK_DEBUG
      std::clog << "GDAL Plugin: featureset, dataset = " << &dataset_ << "\n";
#endif

      query *q = boost::get<query>(&gquery_);
      if(q) {
          return get_feature(*q);
      } else {
          coord2d *p = boost::get<coord2d>(&gquery_);
          if(p) {
              return get_feature_at_point(*p);
          }
      }
      // should never reach here
   }
   return feature_ptr();
}

feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
{
   feature_ptr feature(new Feature(1));

   GDALRasterBand * red = 0;
   GDALRasterBand * green = 0;
   GDALRasterBand * blue = 0;
   GDALRasterBand * alpha = 0;
   GDALRasterBand * grey = 0; 

   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]; // 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);

   // 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::clog << "GDAL Plugin: Raster extent=" << raster_extent << "\n";
   std::clog << "GDAL Plugin: View extent=" << intersect << "\n";
   std::clog << "GDAL Plugin: Query resolution=" << q.resolution() << "\n";
   std::clog << 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 = int(q.resolution() * intersect.width() + 0.5);
      int im_height = int(q.resolution() * intersect.height() + 0.5);

      // case where we need to avoid upsampling so that the
      // image can be later scaled within raster_symbolizer
      if (im_width >= width || im_height >= height)
      {
          im_width = width;
          im_height = height;
      }

      mapnik::ImageData32 image(im_width, im_height);
      image.set(0xffffffff);
         
#ifdef MAPNIK_DEBUG
      std::clog << "GDAL Plugin: Image Size=(" << im_width << "," << im_height << ")\n";
      std::clog << "GDAL Plugin: Reading band " << band_ << "\n";
#endif

      if (band_>0) // we are querying a single band
      {
         float *imageData = (float*)image.getBytes();
         GDALRasterBand * band = dataset_.GetRasterBand(band_);
         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(intersect,image)));
      }
      
      else // working with all bands
      {
         for (int i = 0; i < nbands; ++i)
         {
            GDALRasterBand * band = dataset_.GetRasterBand(i+1);
         
#ifdef MAPNIK_DEBUG
            get_overview_meta(band);  
#endif

            GDALColorInterp color_interp = band->GetColorInterpretation();
            switch (color_interp)
            {
               case GCI_RedBand:
                  red = band;
#ifdef MAPNIK_DEBUG
                  std::clog << "GDAL Plugin: Found red band" << "\n";
#endif
                  break;
               case GCI_GreenBand:
                  green = band;
#ifdef MAPNIK_DEBUG
                  std::clog << "GDAL Plugin: Found green band" << "\n";
#endif
                  break;
               case GCI_BlueBand:
                  blue = band;
#ifdef MAPNIK_DEBUG
                  std::clog << "GDAL Plugin: Found blue band" << "\n";
#endif
                  break;
               case GCI_AlphaBand:
                  alpha = band;
#ifdef MAPNIK_DEBUG
                  std::clog << "GDAL Plugin: Found alpha band" << "\n";
#endif
                  break;
               case GCI_GrayIndex:
                  grey = band;
#ifdef MAPNIK_DEBUG
                  std::clog << "GDAL Plugin: Found gray band" << "\n";
#endif
                  break;
               case GCI_PaletteIndex:
               {
                  grey = band;
#ifdef MAPNIK_DEBUG
                  std::clog << "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::clog << "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::clog << "GDAL Plugin: Color entry RGB (" << ce->c1 << "," <<ce->c2 << "," << ce->c3 << ")\n"; 
#endif
                     }
                  }
                  break;
               }
               case GCI_Undefined:
#ifdef MAPNIK_DEBUG
                  std::clog << "GDAL Plugin: Found undefined band (assumming gray band)" << "\n";
#endif
                  grey = band;
                  break;
               default:
#ifdef MAPNIK_DEBUG
                  std::clog << "GDAL Plugin: band type unknown!" << "\n";
#endif
                  break;
            }
         }

          if (red && green && blue)
         {
#ifdef MAPNIK_DEBUG
            std::clog << "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)
         {
#ifdef MAPNIK_DEBUG
            std::clog << "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)
         {
#ifdef MAPNIK_DEBUG
            std::clog << "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(intersect,image)));
      }
      return feature;
   }
   return feature_ptr();
}


feature_ptr gdal_featureset::get_feature_at_point(mapnik::coord2d const& pt)
{
   if (band_>0) {

      unsigned raster_xsize = dataset_.GetRasterXSize();
      unsigned raster_ysize = dataset_.GetRasterYSize();

      double gt[6];
      dataset_.GetGeoTransform(gt);

      double det = gt[1]*gt[5] - gt[2]*gt[4];
      // subtract half a pixel width & height because gdal coord reference
      // is the top-left corner of a pixel, not the center.
      double X = pt.x - gt[0] - gt[1]/2;
      double Y = pt.y - gt[3] - gt[5]/2;
      double det1 = gt[1]*Y + gt[4]*X;
      double det2 = gt[2]*Y + gt[5]*X;
      unsigned x = det2/det, y = det1/det;

      if(0<=x && x<raster_xsize && 0<=y && y<raster_ysize) {
#ifdef MAPNIK_DEBUG         
         std::clog << boost::format("GDAL Plugin: pt.x=%f pt.y=%f\n") % pt.x % pt.y;
         std::clog << boost::format("GDAL Plugin: x=%f y=%f\n") % x % y;
#endif
         GDALRasterBand * band = dataset_.GetRasterBand(band_);
         int hasNoData;
         double nodata = band->GetNoDataValue(&hasNoData);
         double value;
         band->RasterIO(GF_Read, x, y, 1, 1, &value, 1, 1, GDT_Float64, 0, 0);
         if(!hasNoData || value!=nodata) {
            // construct feature
            feature_ptr feature(new Feature(1));
            geometry2d * point = new point_impl;
            point->move_to(pt.x, pt.y);
            feature->add_geometry(point);
            (*feature)["value"] = value;
            return feature;
         }
      }
   }
   return feature_ptr();
}

void gdal_featureset::get_overview_meta(GDALRasterBand * band)
{
    int band_overviews = band->GetOverviewCount();    
    if (band_overviews > 0)
    {
        std::clog << "GDAL Plugin: "<< band_overviews << " overviews found!" << "\n";
        for (int b = 0; b < band_overviews; b++)
        {
            GDALRasterBand * overview = band->GetOverview (b);
            std::clog << boost::format("GDAL Plugin: Overview=%d Width=%d Height=%d \n")
                                   % b % overview->GetXSize() % overview->GetYSize();
        }
    }
    else
    {
        std::clog << "GDAL Plugin: No overviews found!" << "\n";
    }
    
    int bsx,bsy;
    double scale;
    band->GetBlockSize(&bsx,&bsy);
    scale = band->GetScale();
    std::clog << boost::format("GDAL Plugin: Block=%dx%d Scale=%f Type=%s Color=%s \n") % bsx % bsy % scale
       % GDALGetDataTypeName(band->GetRasterDataType())
       % GDALGetColorInterpretationName(band->GetColorInterpretation()); 
}