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, double filter_factor)
    : dataset_(dataset),
      band_(band),
      gquery_(q),
      filter_factor_(filter_factor),
      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] + (raster_height) * tr[2]; // minx
    double y0 = tr[3] + (raster_height) * tr[5]; // miny
    
    double x1 = tr[0] + (raster_width) * tr[1]; // maxx
    double y1 = tr[3] + (raster_width) * tr[4]; // 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);

    //size of resized output pixel in source image domain
    double margin_x = 1.0/(fabs(dx)*q.resolution());
    double margin_y = 1.0/(fabs(dy)*q.resolution());
    if (margin_x < 1)
        margin_x = 1.0;
    if (margin_y < 1)
        margin_y = 1.0;
    //select minimum raster containing whole box
    int x_off = rint(box.minx() - margin_x);
    int y_off = rint(box.miny() - margin_y);
    int end_x = rint(box.maxx() + margin_x);
    int end_y = rint(box.maxy() + margin_y);
    //clip to available data
    if (x_off < 0)
        x_off = 0;
    if (y_off < 0)
        y_off = 0;
    if (end_x > (int)raster_width)
        end_x = raster_width;
    if (end_y > (int)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 Envelope of returned raster
    Envelope<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";
    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);
	
	// if layer-level filter_factor is set, apply it
	if (filter_factor_)
	{
    	im_width *= filter_factor_;
    	im_height *= filter_factor_;
	}
	// otherwise respect symbolizer level factor applied to query, default of 1.0
	else
	{
      double sym_downsample_factor = q.get_filter_factor();
      im_width *= sym_downsample_factor;
      im_height *= sym_downsample_factor;
	}

	// 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;
	}

	if (im_width > 0 && im_height > 0)
	{
    	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
	typedef std::vector<int,int> pallete;
	
    	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
		// TODO : apply colormap if present
		
    		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()); 
}