/***************************************************************************** * * 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 using mapnik::query; using mapnik::Envelope; using mapnik::Feature; using mapnik::feature_ptr; using mapnik::CoordTransform; gdal_featureset::gdal_featureset(GDALDataset & dataset, query const& q) : dataset_(dataset), query_extent_(q.get_bbox()), first_(true) {} gdal_featureset::~gdal_featureset() {} feature_ptr gdal_featureset::next() { if (first_) { first_ = false; feature_ptr feature(new Feature(1)); GDALRasterBand * red = 0; GDALRasterBand * green = 0; GDALRasterBand * blue = 0; GDALRasterBand * alpha = 0; GDALRasterBand * grey = 0; for (int i = 0; i < dataset_.GetRasterCount() ;++i) { GDALRasterBand * band = dataset_.GetRasterBand(i+1); //int bsx,bsy; //band->GetBlockSize(&bsx,&bsy); //std::cout << boost::format("Block=%dx%d Type=%s Color=%s \n") % bsx % bsy // % GDALGetDataTypeName(band->GetRasterDataType()) // % GDALGetColorInterpretationName(band->GetColorInterpretation()); GDALColorInterp color_interp = band->GetColorInterpretation(); switch (color_interp) { case GCI_RedBand: red = band; break; case GCI_GreenBand: green = band; break; case GCI_BlueBand: blue = band; break; case GCI_AlphaBand: alpha = band; break; case GCI_GrayIndex: grey = band; break; case GCI_Undefined: grey = band; default: //grey = band; break; ; } } unsigned raster_xsize = dataset_.GetRasterXSize(); unsigned raster_ysize = dataset_.GetRasterYSize(); double tr[6]; dataset_.GetGeoTransform(tr); double x0 = tr[0]; double y0 = tr[3]; double x1 = tr[0] + raster_xsize * tr[1] + raster_ysize * tr[2]; double y1 = tr[3] + raster_xsize * tr[4] + raster_ysize * tr[5]; Envelope raster_extent(x0,y0,x1,y1); CoordTransform t (raster_xsize,raster_ysize,raster_extent,0,0); Envelope intersection = raster_extent.intersect(query_extent_); Envelope box = t.forward(intersection); int start_x = int(box.minx()); int start_y = int(box.miny()); int width = int(box.width()); int height = int(box.height()); if (width > 0 && height > 0) { mapnik::ImageData32 image(width,height); image.set(0xffffffff); if (red && green && blue) { red->RasterIO (GF_Read,start_x,start_y,width,height, image.getBytes() + 0, width,height, GDT_Byte,4,4*width); green->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 1, width,height, GDT_Byte,4,4*width); blue->RasterIO (GF_Read,start_x,start_y,width,height, image.getBytes() + 2, width,height, GDT_Byte,4,4*width); if (alpha) { alpha->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 3, width,height, GDT_Byte,4,4*width); } } else if (grey) { grey->RasterIO(GF_Read,start_x,start_y,width,height,image.getBytes() + 0, width, height,GDT_Byte,4,4*width); grey->RasterIO(GF_Read,start_x,start_y,width,height,image.getBytes() + 1, width, height,GDT_Byte,4,4*width); grey->RasterIO(GF_Read,start_x,start_y,width,height,image.getBytes() + 2, width, height,GDT_Byte,4,4*width); if (alpha) { alpha->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 3, width,height, GDT_Byte,4,4*width); } } feature->set_raster(mapnik::raster_ptr(new mapnik::raster(intersection,image))); return feature; } } return feature_ptr(); }