change the method for calculating offsets to pass to GDALRasterIO, as the 'start_x, start_y' implementation added in r1321 was leading to large shifts in 1-banded srtm data without overviews - this new approach works with and without overviews and also with 3 band rbg data. Also added more debug output.
This commit is contained in:
Dane Springmeyer
parent
6d3dcdeac6
commit
27fbc4549d
2 changed files with 132 additions and 77 deletions
|
|
@ -33,6 +33,7 @@ using mapnik::CoordTransform;
|
|||
using mapnik::point_impl;
|
||||
using mapnik::geometry2d;
|
||||
|
||||
#define MAPNIK_DEBUG
|
||||
|
||||
gdal_featureset::gdal_featureset(GDALDataset & dataset, int band, gdal_query q)
|
||||
: dataset_(dataset),
|
||||
|
|
@ -43,7 +44,7 @@ gdal_featureset::gdal_featureset(GDALDataset & dataset, int band, gdal_query q)
|
|||
gdal_featureset::~gdal_featureset()
|
||||
{
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::clog << "closing dataset = " << &dataset_ << "\n";
|
||||
std::clog << "GDAL Plugin: closing dataset = " << &dataset_ << "\n";
|
||||
#endif
|
||||
GDALClose(&dataset_);
|
||||
}
|
||||
|
|
@ -54,7 +55,7 @@ feature_ptr gdal_featureset::next()
|
|||
{
|
||||
first_ = false;
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::clog << "featureset, dataset = " << &dataset_ << "\n";
|
||||
std::clog << "GDAL Plugin: featureset, dataset = " << &dataset_ << "\n";
|
||||
#endif
|
||||
|
||||
query *q = boost::get<query>(&gquery_);
|
||||
|
|
@ -71,7 +72,6 @@ feature_ptr gdal_featureset::next()
|
|||
return feature_ptr();
|
||||
}
|
||||
|
||||
|
||||
feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
|
||||
{
|
||||
feature_ptr feature(new Feature(1));
|
||||
|
|
@ -82,98 +82,85 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
|
|||
GDALRasterBand * alpha = 0;
|
||||
GDALRasterBand * grey = 0;
|
||||
|
||||
unsigned raster_xsize = dataset_.GetRasterXSize();
|
||||
unsigned raster_ysize = dataset_.GetRasterYSize();
|
||||
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];
|
||||
double y0 = tr[3];
|
||||
double x1 = tr[0] + raster_xsize * dx + raster_ysize * tr[2];
|
||||
double y1 = tr[3] + raster_xsize * tr[4] + raster_ysize * dy;
|
||||
Envelope<double> raster_extent(x0,y0,x1,y1);
|
||||
CoordTransform t (raster_xsize,raster_ysize,raster_extent,0,0);
|
||||
Envelope<double> intersection = raster_extent.intersect(q.get_bbox());
|
||||
Envelope<double> box = t.forward(intersection);
|
||||
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);
|
||||
|
||||
int start_x, start_y;
|
||||
if (dx>0)
|
||||
start_x = int(std::floor(intersection.minx()-raster_extent.minx())/dx);
|
||||
else
|
||||
// west-left raster
|
||||
start_x = -int(std::floor(raster_extent.maxx()-intersection.maxx())/dx);
|
||||
if(dy>0)
|
||||
start_y = int(std::floor(intersection.miny()-raster_extent.miny())/dy);
|
||||
else
|
||||
// north-up raster
|
||||
start_y = -int(std::floor(raster_extent.maxy()-intersection.maxy())/dy);
|
||||
|
||||
int width = int(box.width()+0.5);
|
||||
int height = int(box.height()+0.5);
|
||||
// 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::cout << "raster_extent=" << raster_extent << "\n";
|
||||
std::cout << "intersection=" << intersection << "\n";
|
||||
std::cout << boost::format("RasterWidth=%d RasterHeight=%d \n") % raster_xsize % raster_ysize;
|
||||
std::cout << boost::format("StartX=%d StartY=%d Width=%d Height=%d \n") % start_x % start_y % width % height;
|
||||
std::cout << "GDAL Plugin: Raster extent=" << raster_extent << "\n";
|
||||
std::cout << "GDAL Plugin: View extent=" << intersect << "\n";
|
||||
std::cout << "GDAL Plugin: Query resolution=" << q.resolution() << "\n";
|
||||
std::cout << 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 = q.resolution()*intersection.width();
|
||||
int im_height = q.resolution()*intersection.height();
|
||||
int im_width = int(q.resolution() * intersect.width() + 0.5);
|
||||
int im_height = int(q.resolution() * intersect.height() + 0.5);
|
||||
|
||||
mapnik::ImageData32 image(im_width, im_height);
|
||||
image.set(0xffffffff);
|
||||
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "band=" << band_ << "\n";
|
||||
std::cout << "im_size=(" << im_width << "," << im_height << ")\n";
|
||||
std::cout << "box_size=(" << width << "," << height << ")\n";
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "GDAL Plugin: Image Size=(" << im_width << "," << im_height << ")\n";
|
||||
std::cout << "GDAL Plugin: Reading band " << band_ << "\n";
|
||||
#endif
|
||||
|
||||
if(band_>0) // we are querying a single band
|
||||
if (band_>0) // we are querying a single band
|
||||
{
|
||||
float *imageData = (float*)image.getBytes();
|
||||
GDALRasterBand * band = dataset_.GetRasterBand(band_);
|
||||
band->RasterIO(GF_Read, start_x, start_y, width, height,
|
||||
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(intersection,image)));
|
||||
feature->set_raster(mapnik::raster_ptr(new mapnik::raster(intersect,image)));
|
||||
}
|
||||
|
||||
else // working with all bands
|
||||
else // working with all bands
|
||||
{
|
||||
for (int i = 0; i < nbands; ++i)
|
||||
{
|
||||
GDALRasterBand * band = dataset_.GetRasterBand(i+1);
|
||||
|
||||
#ifdef MAPNIK_DEBUG
|
||||
int band_overviews = band->GetOverviewCount();
|
||||
// Note: overviews are picked automatically during GdalDataset::RasterIO below
|
||||
// when nXSize > nPixelSpace and nYSize > nLineSpace
|
||||
|
||||
if (band_overviews > 0)
|
||||
{
|
||||
for (int b = 0; b < band_overviews; b++)
|
||||
{
|
||||
GDALRasterBand * overview = band->GetOverview (b);
|
||||
std::cout << boost::format("Overview=%d Width=%d Height=%d \n")
|
||||
% b % overview->GetXSize() % overview->GetYSize();
|
||||
}
|
||||
}
|
||||
|
||||
int bsx,bsy;
|
||||
double scale;
|
||||
band->GetBlockSize(&bsx,&bsy);
|
||||
scale = band->GetScale();
|
||||
std::cout << boost::format("Block=%dx%d Scale=%f Type=%s Color=%s \n") % bsx % bsy % scale
|
||||
% GDALGetDataTypeName(band->GetRasterDataType())
|
||||
% GDALGetColorInterpretationName(band->GetColorInterpretation());
|
||||
get_overview_meta(band);
|
||||
#endif
|
||||
|
||||
GDALColorInterp color_interp = band->GetColorInterpretation();
|
||||
|
|
@ -181,68 +168,108 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
|
|||
{
|
||||
case GCI_RedBand:
|
||||
red = band;
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "GDAL Plugin: Found red band" << "\n";
|
||||
#endif
|
||||
break;
|
||||
case GCI_GreenBand:
|
||||
green = band;
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "GDAL Plugin: Found green band" << "\n";
|
||||
#endif
|
||||
break;
|
||||
case GCI_BlueBand:
|
||||
blue = band;
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "GDAL Plugin: Found blue band" << "\n";
|
||||
#endif
|
||||
break;
|
||||
case GCI_AlphaBand:
|
||||
alpha = band;
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "GDAL Plugin: Found alpha band" << "\n";
|
||||
#endif
|
||||
break;
|
||||
case GCI_GrayIndex:
|
||||
grey = band;
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "GDAL Plugin: Found gray band" << "\n";
|
||||
#endif
|
||||
break;
|
||||
case GCI_PaletteIndex:
|
||||
{
|
||||
grey = band;
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "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::cout << "Color Table count = " << count << "\n";
|
||||
std::cout << "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::cout << "color entry RGB (" << ce->c1 <<"," <<ce->c2 << "," << ce->c3 << ")\n";
|
||||
std::cout << "GDAL Plugin: Color entry RGB (" << ce->c1 << "," <<ce->c2 << "," << ce->c3 << ")\n";
|
||||
#endif
|
||||
}
|
||||
}
|
||||
break;
|
||||
}
|
||||
case GCI_Undefined:
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "GDAL Plugin: Found undefined band (assumming gray band)" << "\n";
|
||||
#endif
|
||||
grey = band;
|
||||
break;
|
||||
default:
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "GDAL Plugin: band type unknown!" << "\n";
|
||||
#endif
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (red && green && blue)
|
||||
if (red && green && blue)
|
||||
{
|
||||
red->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 0, image.width(),image.height(), GDT_Byte,4,4*image.width());
|
||||
green->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 1, image.width(),image.height(), GDT_Byte,4,4*image.width());
|
||||
blue->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 2, image.width(),image.height(), GDT_Byte,4,4*image.width());
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "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)
|
||||
{
|
||||
grey->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 0, image.width(),image.height(), GDT_Byte,4,4*image.width());
|
||||
grey->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 1, image.width(),image.height(), GDT_Byte,4,4*image.width());
|
||||
grey->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 2, image.width(),image.height(), GDT_Byte,4,4*image.width());
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "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)
|
||||
{
|
||||
alpha->RasterIO(GF_Read,start_x,start_y,width,height, image.getBytes() + 3, image.width(),image.height(), GDT_Byte,4,4*image.width());
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << "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(intersection,image)));
|
||||
feature->set_raster(mapnik::raster_ptr(new mapnik::raster(intersect,image)));
|
||||
}
|
||||
return feature;
|
||||
}
|
||||
|
|
@ -250,7 +277,6 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
|
|||
}
|
||||
|
||||
|
||||
|
||||
feature_ptr gdal_featureset::get_feature_at_point(mapnik::coord2d const& pt)
|
||||
{
|
||||
if (band_>0) {
|
||||
|
|
@ -272,8 +298,8 @@ feature_ptr gdal_featureset::get_feature_at_point(mapnik::coord2d const& pt)
|
|||
|
||||
if(0<=x && x<raster_xsize && 0<=y && y<raster_ysize) {
|
||||
#ifdef MAPNIK_DEBUG
|
||||
std::cout << boost::format("pt.x=%f pt.y=%f\n") % pt.x % pt.y;
|
||||
std::cout << boost::format("x=%f y=%f\n") % x % y;
|
||||
std::cout << boost::format("GDAL Plugin: pt.x=%f pt.y=%f\n") % pt.x % pt.y;
|
||||
std::cout << boost::format("GDAL Plugin: x=%f y=%f\n") % x % y;
|
||||
#endif
|
||||
GDALRasterBand * band = dataset_.GetRasterBand(band_);
|
||||
int hasNoData;
|
||||
|
|
@ -293,3 +319,30 @@ feature_ptr gdal_featureset::get_feature_at_point(mapnik::coord2d const& pt)
|
|||
}
|
||||
return feature_ptr();
|
||||
}
|
||||
|
||||
void gdal_featureset::get_overview_meta(GDALRasterBand * band)
|
||||
{
|
||||
int band_overviews = band->GetOverviewCount();
|
||||
if (band_overviews > 0)
|
||||
{
|
||||
std::cout << "GDAL Plugin: "<< band_overviews << " overviews found!" << "\n";
|
||||
for (int b = 0; b < band_overviews; b++)
|
||||
{
|
||||
GDALRasterBand * overview = band->GetOverview (b);
|
||||
std::cout << boost::format("GDAL Plugin: Overview=%d Width=%d Height=%d \n")
|
||||
% b % overview->GetXSize() % overview->GetYSize();
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cout << "GDAL Plugin: No overviews found!" << "\n";
|
||||
}
|
||||
|
||||
int bsx,bsy;
|
||||
double scale;
|
||||
band->GetBlockSize(&bsx,&bsy);
|
||||
scale = band->GetScale();
|
||||
std::cout << boost::format("GDAL Plugin: Block=%dx%d Scale=%f Type=%s Color=%s \n") % bsx % bsy % scale
|
||||
% GDALGetDataTypeName(band->GetRasterDataType())
|
||||
% GDALGetColorInterpretationName(band->GetColorInterpretation());
|
||||
}
|
||||
|
|
|
|||
|
|
@ -28,6 +28,7 @@
|
|||
#include <boost/variant.hpp>
|
||||
|
||||
class GDALDataset;
|
||||
class GDALRasterBand;
|
||||
|
||||
typedef boost::variant<mapnik::query,mapnik::coord2d> gdal_query;
|
||||
|
||||
|
|
@ -41,6 +42,7 @@ class gdal_featureset : public mapnik::Featureset
|
|||
private:
|
||||
mapnik::feature_ptr get_feature(mapnik::query const& q);
|
||||
mapnik::feature_ptr get_feature_at_point(mapnik::coord2d const& p);
|
||||
void get_overview_meta(GDALRasterBand * band);
|
||||
GDALDataset & dataset_;
|
||||
int band_;
|
||||
gdal_query gquery_;
|
||||
|
|
|
|||
Loading…
Reference in a new issue