/***************************************************************************** * * 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 // boost #include using mapnik::query; using mapnik::coord2d; using mapnik::box2d; 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(&gquery_); if(q) { return get_feature(*q); } else { coord2d *p = boost::get(&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 box2d raster_extent(x0,y0,x1,y1); CoordTransform t(raster_width,raster_height,raster_extent,0,0); box2d intersect = raster_extent.intersect(q.get_bbox()); box2d box = t.forward(intersect); //select minimum raster containing whole box int x_off = static_cast(floor(box.minx())); int y_off = static_cast(floor(box.miny())); int end_x = static_cast(ceil(box.maxx())); int end_y = static_cast(ceil(box.maxy())); //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 box2d of returned raster box2d 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=" << boost::get<0>(q.resolution()) << "," << boost::get<1>(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(boost::get<0>(q.resolution()) * intersect.width() + 0.5); int im_height = int(boost::get<1>(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::image_data_32 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 << "," <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 && xGetNoDataValue(&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()); }