port latest changes from "raster_overzoom_quest" branch.
This commit is contained in:
parent
88fb137061
commit
5d48d98697
2 changed files with 349 additions and 381 deletions
|
@ -187,19 +187,9 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
|
|||
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<double> feature_raster_extent(x_off, y_off, x_off + width, y_off + height);
|
||||
intersect = t.backward(feature_raster_extent);
|
||||
feature_raster_extent = t.backward(feature_raster_extent);
|
||||
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Raster extent=" << raster_extent_;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: View extent=" << intersect;
|
||||
|
@ -208,374 +198,380 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
|
|||
|
||||
if (width > 0 && height > 0)
|
||||
{
|
||||
double width_res = std::get<0>(q.resolution());
|
||||
double height_res = std::get<1>(q.resolution());
|
||||
int im_width = int(width_res * intersect.width() + 0.5);
|
||||
int im_height = int(height_res * intersect.height() + 0.5);
|
||||
|
||||
double filter_factor = q.get_filter_factor();
|
||||
im_width = int(im_width * filter_factor + 0.5);
|
||||
im_height = int(im_height * filter_factor + 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)
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Image Size=(" << width << "," << height << ")";
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Reading band=" << band_;
|
||||
if (band_ > 0) // we are querying a single band
|
||||
{
|
||||
im_width = width;
|
||||
im_height = height;
|
||||
}
|
||||
|
||||
if (im_width > 0 && im_height > 0)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Image Size=(" << im_width << "," << im_height << ")";
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Reading band=" << band_;
|
||||
if (band_ > 0) // we are querying a single band
|
||||
GDALRasterBand * band = dataset_.GetRasterBand(band_);
|
||||
if (band_ > nbands_)
|
||||
{
|
||||
GDALRasterBand * band = dataset_.GetRasterBand(band_);
|
||||
if (band_ > nbands_)
|
||||
std::ostringstream s;
|
||||
s << "GDAL Plugin: " << band_ << " is an invalid band, dataset only has " << nbands_ << "bands";
|
||||
throw datasource_exception(s.str());
|
||||
}
|
||||
GDALDataType band_type = band->GetRasterDataType();
|
||||
switch (band_type)
|
||||
{
|
||||
case GDT_Byte:
|
||||
{
|
||||
mapnik::image_gray8 image(width, height);
|
||||
image.set(std::numeric_limits<std::uint8_t>::max());
|
||||
raster_nodata = band->GetNoDataValue(&raster_has_nodata);
|
||||
raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
image.data(), image.width(), image.height(),
|
||||
GDT_Byte, 0, 0);
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
std::ostringstream s;
|
||||
s << "GDAL Plugin: " << band_ << " is an invalid band, dataset only has " << nbands_ << "bands";
|
||||
throw datasource_exception(s.str());
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
GDALDataType band_type = band->GetRasterDataType();
|
||||
switch (band_type)
|
||||
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(feature_raster_extent, intersect, image, 0.0);
|
||||
// set nodata value to be used in raster colorizer
|
||||
if (nodata_value_) raster->set_nodata(*nodata_value_);
|
||||
else raster->set_nodata(raster_nodata);
|
||||
feature->set_raster(raster);
|
||||
break;
|
||||
}
|
||||
case GDT_Float64:
|
||||
case GDT_Float32:
|
||||
{
|
||||
mapnik::image_gray32f image(width, height);
|
||||
image.set(std::numeric_limits<float>::max());
|
||||
raster_nodata = band->GetNoDataValue(&raster_has_nodata);
|
||||
raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
image.data(), image.width(), image.height(),
|
||||
GDT_Float32, 0, 0);
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
case GDT_Byte:
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(feature_raster_extent, intersect, image, 0.0);
|
||||
// set nodata value to be used in raster colorizer
|
||||
if (nodata_value_) raster->set_nodata(*nodata_value_);
|
||||
else raster->set_nodata(raster_nodata);
|
||||
feature->set_raster(raster);
|
||||
break;
|
||||
}
|
||||
case GDT_UInt16:
|
||||
{
|
||||
mapnik::image_gray16 image(width, height);
|
||||
image.set(std::numeric_limits<std::uint16_t>::max());
|
||||
raster_nodata = band->GetNoDataValue(&raster_has_nodata);
|
||||
raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
image.data(), image.width(), image.height(),
|
||||
GDT_UInt16, 0, 0);
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
mapnik::image_gray8 image(im_width, im_height);
|
||||
image.set(std::numeric_limits<std::uint8_t>::max());
|
||||
raster_nodata = band->GetNoDataValue(&raster_has_nodata);
|
||||
raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
image.data(), image.width(), image.height(),
|
||||
GDT_Byte, 0, 0);
|
||||
if (raster_io_error == CE_Failure)
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(feature_raster_extent, intersect, image, 0.0);
|
||||
// set nodata value to be used in raster colorizer
|
||||
if (nodata_value_) raster->set_nodata(*nodata_value_);
|
||||
else raster->set_nodata(raster_nodata);
|
||||
feature->set_raster(raster);
|
||||
break;
|
||||
}
|
||||
default:
|
||||
case GDT_Int16:
|
||||
{
|
||||
mapnik::image_gray16s image(width, height);
|
||||
image.set(std::numeric_limits<std::int16_t>::max());
|
||||
raster_nodata = band->GetNoDataValue(&raster_has_nodata);
|
||||
raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
image.data(), image.width(), image.height(),
|
||||
GDT_Int16, 0, 0);
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(feature_raster_extent, intersect, image, 0.0);
|
||||
// set nodata value to be used in raster colorizer
|
||||
if (nodata_value_) raster->set_nodata(*nodata_value_);
|
||||
else raster->set_nodata(raster_nodata);
|
||||
feature->set_raster(raster);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
else // working with all bands
|
||||
{
|
||||
mapnik::image_rgba8 image(width, height);
|
||||
image.set(std::numeric_limits<std::uint32_t>::max());
|
||||
for (int i = 0; i < nbands_; ++i)
|
||||
{
|
||||
GDALRasterBand * band = dataset_.GetRasterBand(i + 1);
|
||||
#ifdef MAPNIK_LOG
|
||||
get_overview_meta(band);
|
||||
#endif
|
||||
GDALColorInterp color_interp = band->GetColorInterpretation();
|
||||
switch (color_interp)
|
||||
{
|
||||
case GCI_RedBand:
|
||||
red = band;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found red band";
|
||||
break;
|
||||
case GCI_GreenBand:
|
||||
green = band;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found green band";
|
||||
break;
|
||||
case GCI_BlueBand:
|
||||
blue = band;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found blue band";
|
||||
break;
|
||||
case GCI_AlphaBand:
|
||||
alpha = band;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found alpha band";
|
||||
break;
|
||||
case GCI_GrayIndex:
|
||||
grey = band;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found gray band";
|
||||
break;
|
||||
case GCI_PaletteIndex:
|
||||
{
|
||||
grey = band;
|
||||
#ifdef MAPNIK_LOG
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found gray band, and colortable...";
|
||||
GDALColorTable *color_table = band->GetColorTable();
|
||||
|
||||
if (color_table)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
int count = color_table->GetColorEntryCount();
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Color Table count=" << count;
|
||||
|
||||
for (int j = 0; j < count; j++)
|
||||
{
|
||||
const GDALColorEntry *ce = color_table->GetColorEntry (j);
|
||||
if (! ce) continue;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Color entry RGB=" << ce->c1 << "," <<ce->c2 << "," << ce->c3;
|
||||
}
|
||||
}
|
||||
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(intersect, image, filter_factor);
|
||||
// set nodata value to be used in raster colorizer
|
||||
if (nodata_value_) raster->set_nodata(*nodata_value_);
|
||||
else raster->set_nodata(raster_nodata);
|
||||
feature->set_raster(raster);
|
||||
#endif
|
||||
break;
|
||||
}
|
||||
case GDT_Float64:
|
||||
case GDT_Float32:
|
||||
case GCI_Undefined:
|
||||
#if GDAL_VERSION_NUM <= 1730
|
||||
if (nbands_ == 4)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming alpha band)";
|
||||
alpha = band;
|
||||
}
|
||||
else
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming gray band)";
|
||||
grey = band;
|
||||
}
|
||||
#else
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming gray band)";
|
||||
grey = band;
|
||||
#endif
|
||||
break;
|
||||
default:
|
||||
MAPNIK_LOG_WARN(gdal) << "gdal_featureset: Band type unknown!";
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (red && green && blue)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Processing rgb bands...";
|
||||
raster_nodata = red->GetNoDataValue(&raster_has_nodata);
|
||||
GDALColorTable *color_table = red->GetColorTable();
|
||||
bool has_nodata = nodata_value_ || raster_has_nodata;
|
||||
|
||||
// we can deduce the alpha channel from nodata in the Byte case
|
||||
// by reusing the reading of R,G,B bands directly
|
||||
if (has_nodata && !color_table && red->GetRasterDataType() == GDT_Byte)
|
||||
{
|
||||
mapnik::image_gray32f image(im_width, im_height);
|
||||
image.set(std::numeric_limits<float>::max());
|
||||
raster_nodata = band->GetNoDataValue(&raster_has_nodata);
|
||||
raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
image.data(), image.width(), image.height(),
|
||||
double apply_nodata = nodata_value_ ? *nodata_value_ : raster_nodata;
|
||||
// read the data in and create an alpha channel from the nodata values
|
||||
// TODO - we assume here the nodata value for the red band applies to all bands
|
||||
// more details about this at http://trac.osgeo.org/gdal/ticket/2734
|
||||
float* imageData = (float*)image.bytes();
|
||||
raster_io_error = red->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
imageData, image.width(), image.height(),
|
||||
GDT_Float32, 0, 0);
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
int len = image.width() * image.height();
|
||||
for (int i = 0; i < len; ++i)
|
||||
{
|
||||
if (std::fabs(apply_nodata - imageData[i]) < nodata_tolerance_)
|
||||
{
|
||||
*reinterpret_cast<unsigned *>(&imageData[i]) = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
*reinterpret_cast<unsigned *>(&imageData[i]) = 0xFFFFFFFF;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* Use dataset RasterIO in priority in 99.9% of the cases */
|
||||
if( red->GetBand() == 1 && green->GetBand() == 2 && blue->GetBand() == 3 )
|
||||
{
|
||||
int nBandsToRead = 3;
|
||||
if( alpha != nullptr && alpha->GetBand() == 4 && !raster_has_nodata )
|
||||
{
|
||||
nBandsToRead = 4;
|
||||
alpha = nullptr; // to avoid reading it again afterwards
|
||||
}
|
||||
raster_io_error = dataset_.RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
image.bytes(),
|
||||
image.width(), image.height(), GDT_Byte,
|
||||
nBandsToRead, nullptr,
|
||||
4, 4 * image.width(), 1);
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
raster_io_error = red->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 0,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
raster_io_error = green->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 1,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
raster_io_error = blue->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 2,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
}
|
||||
|
||||
// In the case we skipped initializing the alpha channel
|
||||
if (has_nodata && !color_table && red->GetRasterDataType() == GDT_Byte)
|
||||
{
|
||||
double apply_nodata = nodata_value_ ? *nodata_value_ : raster_nodata;
|
||||
if( apply_nodata >= 0 && apply_nodata <= 255 )
|
||||
{
|
||||
int len = image.width() * image.height();
|
||||
GByte* pabyBytes = (GByte*) image.bytes();
|
||||
for (int i = 0; i < len; ++i)
|
||||
{
|
||||
// TODO - we assume here the nodata value for the red band applies to all bands
|
||||
// more details about this at http://trac.osgeo.org/gdal/ticket/2734
|
||||
if (std::fabs(apply_nodata - pabyBytes[4*i]) < nodata_tolerance_)
|
||||
pabyBytes[4*i + 3] = 0;
|
||||
else
|
||||
pabyBytes[4*i + 3] = 255;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (grey)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Processing gray band...";
|
||||
raster_nodata = grey->GetNoDataValue(&raster_has_nodata);
|
||||
GDALColorTable* color_table = grey->GetColorTable();
|
||||
bool has_nodata = nodata_value_ || raster_has_nodata;
|
||||
if (!color_table && has_nodata)
|
||||
{
|
||||
double apply_nodata = nodata_value_ ? *nodata_value_ : raster_nodata;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: applying nodata value for layer=" << apply_nodata;
|
||||
// first read the data in and create an alpha channel from the nodata values
|
||||
float* imageData = (float*)image.bytes();
|
||||
raster_io_error = grey->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
imageData, image.width(), image.height(),
|
||||
GDT_Float32, 0, 0);
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(intersect, image, filter_factor);
|
||||
// set nodata value to be used in raster colorizer
|
||||
if (nodata_value_) raster->set_nodata(*nodata_value_);
|
||||
else raster->set_nodata(raster_nodata);
|
||||
feature->set_raster(raster);
|
||||
break;
|
||||
}
|
||||
case GDT_UInt16:
|
||||
{
|
||||
mapnik::image_gray16 image(im_width, im_height);
|
||||
image.set(std::numeric_limits<std::uint16_t>::max());
|
||||
raster_nodata = band->GetNoDataValue(&raster_has_nodata);
|
||||
raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
image.data(), image.width(), image.height(),
|
||||
GDT_UInt16, 0, 0);
|
||||
if (raster_io_error == CE_Failure)
|
||||
int len = image.width() * image.height();
|
||||
for (int i = 0; i < len; ++i)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(intersect, image, filter_factor);
|
||||
// set nodata value to be used in raster colorizer
|
||||
if (nodata_value_) raster->set_nodata(*nodata_value_);
|
||||
else raster->set_nodata(raster_nodata);
|
||||
feature->set_raster(raster);
|
||||
break;
|
||||
}
|
||||
default:
|
||||
case GDT_Int16:
|
||||
{
|
||||
mapnik::image_gray16s image(im_width, im_height);
|
||||
image.set(std::numeric_limits<std::int16_t>::max());
|
||||
raster_nodata = band->GetNoDataValue(&raster_has_nodata);
|
||||
raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
image.data(), image.width(), image.height(),
|
||||
GDT_Int16, 0, 0);
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(intersect, image, filter_factor);
|
||||
// set nodata value to be used in raster colorizer
|
||||
if (nodata_value_) raster->set_nodata(*nodata_value_);
|
||||
else raster->set_nodata(raster_nodata);
|
||||
feature->set_raster(raster);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
else // working with all bands
|
||||
{
|
||||
mapnik::image_rgba8 image(im_width, im_height);
|
||||
image.set(std::numeric_limits<std::uint32_t>::max());
|
||||
for (int i = 0; i < nbands_; ++i)
|
||||
{
|
||||
GDALRasterBand * band = dataset_.GetRasterBand(i + 1);
|
||||
#ifdef MAPNIK_LOG
|
||||
get_overview_meta(band);
|
||||
#endif
|
||||
GDALColorInterp color_interp = band->GetColorInterpretation();
|
||||
switch (color_interp)
|
||||
{
|
||||
case GCI_RedBand:
|
||||
red = band;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found red band";
|
||||
break;
|
||||
case GCI_GreenBand:
|
||||
green = band;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found green band";
|
||||
break;
|
||||
case GCI_BlueBand:
|
||||
blue = band;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found blue band";
|
||||
break;
|
||||
case GCI_AlphaBand:
|
||||
alpha = band;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found alpha band";
|
||||
break;
|
||||
case GCI_GrayIndex:
|
||||
grey = band;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found gray band";
|
||||
break;
|
||||
case GCI_PaletteIndex:
|
||||
{
|
||||
grey = band;
|
||||
#ifdef MAPNIK_LOG
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found gray band, and colortable...";
|
||||
GDALColorTable *color_table = band->GetColorTable();
|
||||
|
||||
if (color_table)
|
||||
if (std::fabs(apply_nodata - imageData[i]) < nodata_tolerance_)
|
||||
{
|
||||
int count = color_table->GetColorEntryCount();
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Color Table count=" << count;
|
||||
|
||||
for (int j = 0; j < count; j++)
|
||||
{
|
||||
const GDALColorEntry *ce = color_table->GetColorEntry (j);
|
||||
if (! ce) continue;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Color entry RGB=" << ce->c1 << "," <<ce->c2 << "," << ce->c3;
|
||||
}
|
||||
}
|
||||
#endif
|
||||
break;
|
||||
}
|
||||
case GCI_Undefined:
|
||||
#if GDAL_VERSION_NUM <= 1730
|
||||
if (nbands_ == 4)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming alpha band)";
|
||||
alpha = band;
|
||||
*reinterpret_cast<unsigned *>(&imageData[i]) = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming gray band)";
|
||||
grey = band;
|
||||
*reinterpret_cast<unsigned *>(&imageData[i]) = 0xFFFFFFFF;
|
||||
}
|
||||
#else
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming gray band)";
|
||||
grey = band;
|
||||
#endif
|
||||
break;
|
||||
default:
|
||||
MAPNIK_LOG_WARN(gdal) << "gdal_featureset: Band type unknown!";
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (red && green && blue)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Processing rgb bands...";
|
||||
raster_nodata = red->GetNoDataValue(&raster_has_nodata);
|
||||
GDALColorTable *color_table = red->GetColorTable();
|
||||
bool has_nodata = nodata_value_ || raster_has_nodata;
|
||||
|
||||
// we can deduce the alpha channel from nodata in the Byte case
|
||||
// by reusing the reading of R,G,B bands directly
|
||||
if (has_nodata && !color_table && red->GetRasterDataType() == GDT_Byte)
|
||||
raster_io_error = grey->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 0,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
|
||||
raster_io_error = grey->RasterIO(GF_Read,x_off, y_off, width, height, image.bytes() + 1,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
|
||||
raster_io_error = grey->RasterIO(GF_Read,x_off, y_off, width, height, image.bytes() + 2,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
|
||||
if (color_table)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Loading color table...";
|
||||
for (unsigned y = 0; y < image.height(); ++y)
|
||||
{
|
||||
double apply_nodata = nodata_value_ ? *nodata_value_ : raster_nodata;
|
||||
// read the data in and create an alpha channel from the nodata values
|
||||
// TODO - we assume here the nodata value for the red band applies to all bands
|
||||
// more details about this at http://trac.osgeo.org/gdal/ticket/2734
|
||||
float* imageData = (float*)image.bytes();
|
||||
raster_io_error = red->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
imageData, image.width(), image.height(),
|
||||
GDT_Float32, 0, 0);
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
int len = image.width() * image.height();
|
||||
for (int i = 0; i < len; ++i)
|
||||
unsigned int* row = image.get_row(y);
|
||||
for (unsigned x = 0; x < image.width(); ++x)
|
||||
{
|
||||
if (std::fabs(apply_nodata - imageData[i]) < nodata_tolerance_)
|
||||
unsigned value = row[x] & 0xff;
|
||||
const GDALColorEntry *ce = color_table->GetColorEntry(value);
|
||||
if (ce)
|
||||
{
|
||||
*reinterpret_cast<unsigned *>(&imageData[i]) = 0;
|
||||
row[x] = (ce->c4 << 24)| (ce->c3 << 16) | (ce->c2 << 8) | (ce->c1) ;
|
||||
}
|
||||
else
|
||||
{
|
||||
*reinterpret_cast<unsigned *>(&imageData[i]) = 0xFFFFFFFF;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* Use dataset RasterIO in priority in 99.9% of the cases */
|
||||
if( red->GetBand() == 1 && green->GetBand() == 2 && blue->GetBand() == 3 )
|
||||
{
|
||||
int nBandsToRead = 3;
|
||||
if( alpha != nullptr && alpha->GetBand() == 4 && !raster_has_nodata )
|
||||
{
|
||||
nBandsToRead = 4;
|
||||
alpha = nullptr; // to avoid reading it again afterwards
|
||||
}
|
||||
raster_io_error = dataset_.RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
image.bytes(),
|
||||
image.width(), image.height(), GDT_Byte,
|
||||
nBandsToRead, nullptr,
|
||||
4, 4 * image.width(), 1);
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
raster_io_error = red->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 0,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
raster_io_error = green->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 1,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
raster_io_error = blue->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 2,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
}
|
||||
|
||||
// In the case we skipped initializing the alpha channel
|
||||
if (has_nodata && !color_table && red->GetRasterDataType() == GDT_Byte)
|
||||
{
|
||||
double apply_nodata = nodata_value_ ? *nodata_value_ : raster_nodata;
|
||||
if( apply_nodata >= 0 && apply_nodata <= 255 )
|
||||
{
|
||||
int len = image.width() * image.height();
|
||||
GByte* pabyBytes = (GByte*) image.bytes();
|
||||
for (int i = 0; i < len; ++i)
|
||||
{
|
||||
// TODO - we assume here the nodata value for the red band applies to all bands
|
||||
// more details about this at http://trac.osgeo.org/gdal/ticket/2734
|
||||
if (std::fabs(apply_nodata - pabyBytes[4*i]) < nodata_tolerance_)
|
||||
pabyBytes[4*i + 3] = 0;
|
||||
else
|
||||
pabyBytes[4*i + 3] = 255;
|
||||
// make lacking color entry fully alpha
|
||||
// note - gdal_translate makes black
|
||||
row[x] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (grey)
|
||||
}
|
||||
if (alpha)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: processing alpha band...";
|
||||
if (!raster_has_nodata)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Processing gray band...";
|
||||
raster_nodata = grey->GetNoDataValue(&raster_has_nodata);
|
||||
GDALColorTable* color_table = grey->GetColorTable();
|
||||
bool has_nodata = nodata_value_ || raster_has_nodata;
|
||||
if (!color_table && has_nodata)
|
||||
{
|
||||
double apply_nodata = nodata_value_ ? *nodata_value_ : raster_nodata;
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: applying nodata value for layer=" << apply_nodata;
|
||||
// first read the data in and create an alpha channel from the nodata values
|
||||
float* imageData = (float*)image.bytes();
|
||||
raster_io_error = grey->RasterIO(GF_Read, x_off, y_off, width, height,
|
||||
imageData, image.width(), image.height(),
|
||||
GDT_Float32, 0, 0);
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
int len = image.width() * image.height();
|
||||
for (int i = 0; i < len; ++i)
|
||||
{
|
||||
if (std::fabs(apply_nodata - imageData[i]) < nodata_tolerance_)
|
||||
{
|
||||
*reinterpret_cast<unsigned *>(&imageData[i]) = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
*reinterpret_cast<unsigned *>(&imageData[i]) = 0xFFFFFFFF;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
raster_io_error = grey->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 0,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
raster_io_error = alpha->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 3,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
|
||||
raster_io_error = grey->RasterIO(GF_Read,x_off, y_off, width, height, image.bytes() + 1,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
|
||||
raster_io_error = grey->RasterIO(GF_Read,x_off, y_off, width, height, image.bytes() + 2,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
|
||||
if (raster_io_error == CE_Failure)
|
||||
{
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
|
||||
if (color_table)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Loading color table...";
|
||||
for (unsigned y = 0; y < image.height(); ++y)
|
||||
{
|
||||
unsigned int* row = image.get_row(y);
|
||||
for (unsigned x = 0; x < image.width(); ++x)
|
||||
{
|
||||
unsigned value = row[x] & 0xff;
|
||||
const GDALColorEntry *ce = color_table->GetColorEntry(value);
|
||||
if (ce)
|
||||
{
|
||||
row[x] = (ce->c4 << 24)| (ce->c3 << 16) | (ce->c2 << 8) | (ce->c1) ;
|
||||
}
|
||||
else
|
||||
{
|
||||
// make lacking color entry fully alpha
|
||||
// note - gdal_translate makes black
|
||||
row[x] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if (alpha)
|
||||
else
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: processing alpha band...";
|
||||
MAPNIK_LOG_WARN(gdal) << "warning: nodata value (" << raster_nodata << ") used to set transparency instead of alpha band";
|
||||
}
|
||||
}
|
||||
else if( dataset_.GetRasterCount() > 0 && dataset_.GetRasterBand(1) )
|
||||
{
|
||||
// Check if we have a non-alpha mask band (for example a TIFF internal mask)
|
||||
int flags = dataset_.GetRasterBand(1)->GetMaskFlags();
|
||||
GDALRasterBand* mask = 0;
|
||||
if (flags == GMF_PER_DATASET)
|
||||
{
|
||||
mask = dataset_.GetRasterBand(1)->GetMaskBand();
|
||||
}
|
||||
if (mask)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: found and processing mask band...";
|
||||
if (!raster_has_nodata)
|
||||
{
|
||||
raster_io_error = alpha->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 3,
|
||||
raster_io_error = mask->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 3,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
|
@ -583,48 +579,22 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
|
|||
}
|
||||
else
|
||||
{
|
||||
MAPNIK_LOG_WARN(gdal) << "warning: nodata value (" << raster_nodata << ") used to set transparency instead of alpha band";
|
||||
MAPNIK_LOG_WARN(gdal) << "warning: nodata value (" << raster_nodata << ") used to set transparency instead of mask band";
|
||||
}
|
||||
}
|
||||
else if( dataset_.GetRasterCount() > 0 && dataset_.GetRasterBand(1) )
|
||||
{
|
||||
// Check if we have a non-alpha mask band (for example a TIFF internal mask)
|
||||
int flags = dataset_.GetRasterBand(1)->GetMaskFlags();
|
||||
GDALRasterBand* mask = 0;
|
||||
if (flags == GMF_PER_DATASET)
|
||||
{
|
||||
mask = dataset_.GetRasterBand(1)->GetMaskBand();
|
||||
}
|
||||
if (mask)
|
||||
{
|
||||
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: found and processing mask band...";
|
||||
if (!raster_has_nodata)
|
||||
{
|
||||
raster_io_error = mask->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 3,
|
||||
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
|
||||
if (raster_io_error == CE_Failure) {
|
||||
throw datasource_exception(CPLGetLastErrorMsg());
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
MAPNIK_LOG_WARN(gdal) << "warning: nodata value (" << raster_nodata << ") used to set transparency instead of mask band";
|
||||
}
|
||||
}
|
||||
}
|
||||
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(intersect, image, filter_factor);
|
||||
// set nodata value to be used in raster colorizer
|
||||
if (nodata_value_) raster->set_nodata(*nodata_value_);
|
||||
else raster->set_nodata(raster_nodata);
|
||||
feature->set_raster(raster);
|
||||
}
|
||||
// report actual/original source nodata in feature attributes
|
||||
if (raster_has_nodata)
|
||||
{
|
||||
feature->put("nodata",raster_nodata);
|
||||
}
|
||||
return feature;
|
||||
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(feature_raster_extent, intersect, image, 0.0);
|
||||
// set nodata value to be used in raster colorizer
|
||||
if (nodata_value_) raster->set_nodata(*nodata_value_);
|
||||
else raster->set_nodata(raster_nodata);
|
||||
feature->set_raster(raster);
|
||||
}
|
||||
// report actual/original source nodata in feature attributes
|
||||
if (raster_has_nodata)
|
||||
{
|
||||
feature->put("nodata",raster_nodata);
|
||||
}
|
||||
return feature;
|
||||
}
|
||||
return feature_ptr();
|
||||
}
|
||||
|
|
|
@ -96,6 +96,8 @@ feature_ptr raster_featureset<LookupPolicy>::next()
|
|||
int end_y = static_cast<int>(std::ceil(ext.maxy()));
|
||||
|
||||
// clip to available data
|
||||
if (x_off >= image_width) x_off = image_width - 1;
|
||||
if (y_off >= image_width) y_off = image_width - 1;
|
||||
if (x_off < 0) x_off = 0;
|
||||
if (y_off < 0) y_off = 0;
|
||||
if (end_x > image_width) end_x = image_width;
|
||||
|
@ -103,12 +105,8 @@ feature_ptr raster_featureset<LookupPolicy>::next()
|
|||
|
||||
int width = end_x - x_off;
|
||||
int height = end_y - y_off;
|
||||
if (width < 1) {
|
||||
width = 1;
|
||||
}
|
||||
if (height < 1) {
|
||||
height = 1;
|
||||
}
|
||||
if (width < 1) width = 1;
|
||||
if (height < 1) height = 1;
|
||||
|
||||
// calculate actual box2d of returned raster
|
||||
box2d<double> feature_raster_extent(rem.minx() + x_off,
|
||||
|
|
Loading…
Reference in a new issue