sij/mapnik
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge pull request #3668 from mapnik/raster_overzoom_quest

Browse source 
Raster Overzoom Quest
This commit is contained in:
Blake Thompson 2017-05-12 11:48:34 -05:00 committed by GitHub
commit 5911fe6374
8 changed files with 421 additions and 425 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
bootstrap.sh
View file

@ -11,7 +11,7 @@ todo
- shrink icu data - shrink icu data
' '
MASON_VERSION="3c6df04" MASON_VERSION="v0.10.0"
function setup_mason() { function setup_mason() {
if [[ ! -d ./.mason ]]; then if [[ ! -d ./.mason ]]; then

1
demo/viewer/viewer.pro
View file

@ -1,7 +1,6 @@
###################################################################### ######################################################################
# Mapnik viewer - Copyright (C) 2007 Artem Pavlenko # Mapnik viewer - Copyright (C) 2007 Artem Pavlenko
###################################################################### ######################################################################
QMAKE_MAC_SDK = macosx10.11
TEMPLATE = app TEMPLATE = app
QT += core gui widgets QT += core gui widgets
QMAKE_CXX = $$system(mapnik-config --cxx) QMAKE_CXX = $$system(mapnik-config --cxx)

13
include/mapnik/raster.hpp
View file

@ -40,15 +40,27 @@ class raster : private util::noncopyable
{ {
public: public:
box2d<double> ext_; box2d<double> ext_;
box2d<double> query_ext_;
image_any data_; image_any data_;
double filter_factor_; double filter_factor_;
boost::optional<double> nodata_; boost::optional<double> nodata_;
template <typename ImageData>
raster(box2d<double> const& ext,
box2d<double> const& query_ext,
ImageData && data,
double filter_factor)
: ext_(ext),
query_ext_(query_ext),
data_(std::move(data)),
filter_factor_(filter_factor) {}
template <typename ImageData> template <typename ImageData>
raster(box2d<double> const& ext, raster(box2d<double> const& ext,
ImageData && data, ImageData && data,
double filter_factor) double filter_factor)
: ext_(ext), : ext_(ext),
query_ext_(ext),
data_(std::move(data)), data_(std::move(data)),
filter_factor_(filter_factor) {} filter_factor_(filter_factor) {}
@ -71,7 +83,6 @@ public:
{ {
filter_factor_ = factor; filter_factor_ = factor;
} }
}; };
} }

49
include/mapnik/renderer_common/process_raster_symbolizer.hpp
View file

@ -53,6 +53,7 @@ struct image_dispatcher
image_dispatcher(int start_x, int start_y, image_dispatcher(int start_x, int start_y,
int width, int height, int width, int height,
double scale_x, double scale_y, double scale_x, double scale_y,
double offset_x, double offset_y,
scaling_method_e method, double filter_factor, scaling_method_e method, double filter_factor,
double opacity, composite_mode_e comp_op, double opacity, composite_mode_e comp_op,
raster_symbolizer const& sym, feature_impl const& feature, raster_symbolizer const& sym, feature_impl const& feature,
@ -63,15 +64,17 @@ struct image_dispatcher
height_(height), height_(height),
scale_x_(scale_x), scale_x_(scale_x),
scale_y_(scale_y), scale_y_(scale_y),
offset_x_(offset_x),
offset_y_(offset_y),
method_(method), method_(method),
filter_factor_(filter_factor), filter_factor_(filter_factor),
opacity_(opacity), opacity_(opacity),
comp_op_(comp_op), comp_op_(comp_op),
sym_(sym), sym_(sym),
feature_(feature), feature_(feature),
composite_(composite), composite_(composite),
nodata_(nodata), nodata_(nodata),
need_scaling_(need_scaling) {} need_scaling_(need_scaling) {}
void operator() (image_null const&) const {} //no-op void operator() (image_null const&) const {} //no-op
void operator() (image_rgba8 const& data_in) const void operator() (image_rgba8 const& data_in) const
@ -79,7 +82,7 @@ struct image_dispatcher
if (need_scaling_) if (need_scaling_)
{ {
image_rgba8 data_out(width_, height_, true, true); image_rgba8 data_out(width_, height_, true, true);
scale_image_agg(data_out, data_in, method_, scale_x_, scale_y_, 0.0, 0.0, filter_factor_, nodata_); scale_image_agg(data_out, data_in, method_, scale_x_, scale_y_, offset_x_, offset_y_, filter_factor_, nodata_);
composite_(data_out, comp_op_, opacity_, start_x_, start_y_); composite_(data_out, comp_op_, opacity_, start_x_, start_y_);
} }
else else
@ -97,7 +100,7 @@ struct image_dispatcher
if (need_scaling_) if (need_scaling_)
{ {
image_type data_out(width_, height_); image_type data_out(width_, height_);
scale_image_agg(data_out, data_in, method_, scale_x_, scale_y_, 0.0, 0.0, filter_factor_, nodata_); scale_image_agg(data_out, data_in, method_, scale_x_, scale_y_, offset_x_, offset_y_, filter_factor_, nodata_);
if (colorizer) colorizer->colorize(dst, data_out, nodata_, feature_); if (colorizer) colorizer->colorize(dst, data_out, nodata_, feature_);
} }
else else
@ -114,6 +117,8 @@ private:
int height_; int height_;
double scale_x_; double scale_x_;
double scale_y_; double scale_y_;
double offset_x_;
double offset_y_;
scaling_method_e method_; scaling_method_e method_;
double filter_factor_; double filter_factor_;
double opacity_; double opacity_;
@ -210,12 +215,18 @@ void render_raster_symbolizer(raster_symbolizer const& sym,
if (source) if (source)
{ {
box2d<double> target_ext = box2d<double>(source->ext_); box2d<double> target_ext = box2d<double>(source->ext_);
prj_trans.backward(target_ext, PROJ_ENVELOPE_POINTS); box2d<double> target_query_ext = box2d<double>(source->query_ext_);
if (!prj_trans.equal())
{
prj_trans.backward(target_ext, PROJ_ENVELOPE_POINTS);
prj_trans.backward(target_query_ext, PROJ_ENVELOPE_POINTS);
}
box2d<double> ext = common.t_.forward(target_ext); box2d<double> ext = common.t_.forward(target_ext);
int start_x = static_cast<int>(std::floor(ext.minx()+.5)); box2d<double> query_ext = common.t_.forward(target_query_ext);
int start_y = static_cast<int>(std::floor(ext.miny()+.5)); int start_x = static_cast<int>(std::floor(query_ext.minx()+.5));
int end_x = static_cast<int>(std::floor(ext.maxx()+.5)); int start_y = static_cast<int>(std::floor(query_ext.miny()+.5));
int end_y = static_cast<int>(std::floor(ext.maxy()+.5)); int end_x = static_cast<int>(std::floor(query_ext.maxx()+.5));
int end_y = static_cast<int>(std::floor(query_ext.maxy()+.5));
int raster_width = end_x - start_x; int raster_width = end_x - start_x;
int raster_height = end_y - start_y; int raster_height = end_y - start_y;
if (raster_width > 0 && raster_height > 0) if (raster_width > 0 && raster_height > 0)
@ -236,17 +247,20 @@ void render_raster_symbolizer(raster_symbolizer const& sym,
if (!prj_trans.equal()) if (!prj_trans.equal())
{ {
double offset_x = ext.minx() - start_x; // This path does not currently work and is still being figured out.
double offset_y = ext.miny() - start_y; double offset_x = query_ext.minx() - start_x;
double offset_y = query_ext.miny() - start_y;
unsigned mesh_size = static_cast<unsigned>(get<value_integer>(sym,keys::mesh_size,feature, common.vars_, 16)); unsigned mesh_size = static_cast<unsigned>(get<value_integer>(sym,keys::mesh_size,feature, common.vars_, 16));
detail::image_warp_dispatcher<F> dispatcher(prj_trans, start_x, start_y, raster_width, raster_height, detail::image_warp_dispatcher<F> dispatcher(prj_trans, start_x, start_y, raster_width, raster_height,
target_ext, source->ext_, offset_x, offset_y, mesh_size, target_query_ext, source->ext_, offset_x, offset_y, mesh_size,
scaling_method, source->get_filter_factor(), scaling_method, source->get_filter_factor(),
opacity, comp_op, sym, feature, composite, source->nodata()); opacity, comp_op, sym, feature, composite, source->nodata());
util::apply_visitor(dispatcher, source->data_); util::apply_visitor(dispatcher, source->data_);
} }
else else
{ {
double offset_x = query_ext.minx() - ext.minx();
double offset_y = query_ext.miny() - ext.miny();
double image_ratio_x = ext.width() / source->data_.width(); double image_ratio_x = ext.width() / source->data_.width();
double image_ratio_y = ext.height() / source->data_.height(); double image_ratio_y = ext.height() / source->data_.height();
double eps = 1e-5; double eps = 1e-5;
@ -256,6 +270,7 @@ void render_raster_symbolizer(raster_symbolizer const& sym,
(std::abs(start_y) > eps); (std::abs(start_y) > eps);
detail::image_dispatcher<F> dispatcher(start_x, start_y, raster_width, raster_height, detail::image_dispatcher<F> dispatcher(start_x, start_y, raster_width, raster_height,
image_ratio_x, image_ratio_y, image_ratio_x, image_ratio_y,
offset_x, offset_y,
scaling_method, source->get_filter_factor(), scaling_method, source->get_filter_factor(),
opacity, comp_op, sym, feature, composite, source->nodata(), scale); opacity, comp_op, sym, feature, composite, source->nodata(), scale);
util::apply_visitor(dispatcher, source->data_); util::apply_visitor(dispatcher, source->data_);

720
plugins/input/gdal/gdal_featureset.cpp
View file

@ -187,19 +187,9 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
int width = end_x - x_off; int width = end_x - x_off;
int height = end_y - y_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 //calculate actual box2d of returned raster
box2d<double> feature_raster_extent(x_off, y_off, x_off + width, y_off + height); 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: Raster extent=" << raster_extent_;
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: View extent=" << intersect; 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) if (width > 0 && height > 0)
{ {
double width_res = std::get<0>(q.resolution()); MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Image Size=(" << width << "," << height << ")";
double height_res = std::get<1>(q.resolution()); MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Reading band=" << band_;
int im_width = int(width_res * intersect.width() + 0.5); if (band_ > 0) // we are querying a single band
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)
{ {
im_width = width; GDALRasterBand * band = dataset_.GetRasterBand(band_);
im_height = height; if (band_ > nbands_)
}
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_); std::ostringstream s;
if (band_ > nbands_) 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; throw datasource_exception(CPLGetLastErrorMsg());
s << "GDAL Plugin: " << band_ << " is an invalid band, dataset only has " << nbands_ << "bands";
throw datasource_exception(s.str());
} }
GDALDataType band_type = band->GetRasterDataType(); mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(feature_raster_extent, intersect, image, 0.0);
switch (band_type) // 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); throw datasource_exception(CPLGetLastErrorMsg());
image.set(std::numeric_limits<std::uint8_t>::max()); }
raster_nodata = band->GetNoDataValue(&raster_has_nodata); mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(feature_raster_extent, intersect, image, 0.0);
raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height, // set nodata value to be used in raster colorizer
image.data(), image.width(), image.height(), if (nodata_value_) raster->set_nodata(*nodata_value_);
GDT_Byte, 0, 0); else raster->set_nodata(raster_nodata);
if (raster_io_error == CE_Failure) 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); #endif
// 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; break;
} }
case GDT_Float64: case GCI_Undefined:
case GDT_Float32: #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); double apply_nodata = nodata_value_ ? *nodata_value_ : raster_nodata;
image.set(std::numeric_limits<float>::max()); // read the data in and create an alpha channel from the nodata values
raster_nodata = band->GetNoDataValue(&raster_has_nodata); // TODO - we assume here the nodata value for the red band applies to all bands
raster_io_error = band->RasterIO(GF_Read, x_off, y_off, width, height, // more details about this at http://trac.osgeo.org/gdal/ticket/2734
image.data(), image.width(), image.height(), 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); GDT_Float32, 0, 0);
if (raster_io_error == CE_Failure) if (raster_io_error == CE_Failure)
{ {
throw datasource_exception(CPLGetLastErrorMsg()); throw datasource_exception(CPLGetLastErrorMsg());
} }
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(intersect, image, filter_factor); int len = image.width() * image.height();
// set nodata value to be used in raster colorizer for (int i = 0; i < len; ++i)
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)
{ {
throw datasource_exception(CPLGetLastErrorMsg()); if (std::fabs(apply_nodata - imageData[i]) < nodata_tolerance_)
}
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)
{ {
int count = color_table->GetColorEntryCount(); *reinterpret_cast<unsigned *>(&imageData[i]) = 0;
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;
} }
else else
{ {
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming gray band)"; *reinterpret_cast<unsigned *>(&imageData[i]) = 0xFFFFFFFF;
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 raster_io_error = grey->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 0,
// by reusing the reading of R,G,B bands directly image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
if (has_nodata && !color_table && red->GetRasterDataType() == GDT_Byte) 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; unsigned int* row = image.get_row(y);
// read the data in and create an alpha channel from the nodata values for (unsigned x = 0; x < image.width(); ++x)
// 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_) 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 else
{ {
*reinterpret_cast<unsigned *>(&imageData[i]) = 0xFFFFFFFF; // make lacking color entry fully alpha
} // note - gdal_translate makes black
} row[x] = 0;
}
/* 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) }
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_io_error = alpha->RasterIO(GF_Read, x_off, y_off, width, height, image.bytes() + 3,
raster_nodata = grey->GetNoDataValue(&raster_has_nodata); image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
GDALColorTable* color_table = grey->GetColorTable(); if (raster_io_error == CE_Failure) {
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)
{
throw datasource_exception(CPLGetLastErrorMsg()); 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) 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()); image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
if (raster_io_error == CE_Failure) { if (raster_io_error == CE_Failure) {
throw datasource_exception(CPLGetLastErrorMsg()); throw datasource_exception(CPLGetLastErrorMsg());
@ -583,48 +579,22 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
} }
else 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 mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(feature_raster_extent, intersect, image, 0.0);
if (raster_has_nodata) // set nodata value to be used in raster colorizer
{ if (nodata_value_) raster->set_nodata(*nodata_value_);
feature->put("nodata",raster_nodata); else raster->set_nodata(raster_nodata);
} feature->set_raster(raster);
return feature;
} }
// report actual/original source nodata in feature attributes
if (raster_has_nodata)
{
feature->put("nodata",raster_nodata);
}
return feature;
} }
return feature_ptr(); return feature_ptr();
} }

49
plugins/input/raster/raster_featureset.cpp
View file

@ -89,33 +89,34 @@ feature_ptr raster_featureset<LookupPolicy>::next()
box2d<double> intersect = bbox_.intersect(curIter_->envelope()); box2d<double> intersect = bbox_.intersect(curIter_->envelope());
box2d<double> ext = t.forward(intersect); box2d<double> ext = t.forward(intersect);
box2d<double> rem = policy_.transform(ext); box2d<double> rem = policy_.transform(ext);
if (ext.width() > 0.5 && ext.height() > 0.5 ) // select minimum raster containing whole ext
{ int x_off = static_cast<int>(std::floor(ext.minx()));
// select minimum raster containing whole ext int y_off = static_cast<int>(std::floor(ext.miny()));
int x_off = static_cast<int>(std::floor(ext.minx())); int end_x = static_cast<int>(std::ceil(ext.maxx()));
int y_off = static_cast<int>(std::floor(ext.miny())); int end_y = static_cast<int>(std::ceil(ext.maxy()));
int end_x = static_cast<int>(std::ceil(ext.maxx()));
int end_y = static_cast<int>(std::ceil(ext.maxy()));
// clip to available data // clip to available data
if (x_off < 0) x_off = 0; if (x_off >= image_width) x_off = image_width - 1;
if (y_off < 0) y_off = 0; if (y_off >= image_width) y_off = image_width - 1;
if (end_x > image_width) end_x = image_width; if (x_off < 0) x_off = 0;
if (end_y > image_height) end_y = image_height; if (y_off < 0) y_off = 0;
if (end_x > image_width) end_x = image_width;
if (end_y > image_height) end_y = image_height;
int width = end_x - x_off; int width = end_x - x_off;
int height = end_y - y_off; int height = end_y - y_off;
if (width < 1) width = 1;
if (height < 1) height = 1;
// calculate actual box2d of returned raster // calculate actual box2d of returned raster
box2d<double> feature_raster_extent(rem.minx() + x_off, box2d<double> feature_raster_extent(rem.minx() + x_off,
rem.miny() + y_off, rem.miny() + y_off,
rem.maxx() + x_off + width, rem.maxx() + x_off + width,
rem.maxy() + y_off + height); rem.maxy() + y_off + height);
intersect = t.backward(feature_raster_extent); feature_raster_extent = t.backward(feature_raster_extent);
mapnik::image_any data = reader->read(x_off, y_off, width, height); mapnik::image_any data = reader->read(x_off, y_off, width, height);
mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(intersect, std::move(data), 1.0); mapnik::raster_ptr raster = std::make_shared<mapnik::raster>(feature_raster_extent, intersect, std::move(data), 1.0);
feature->set_raster(raster); feature->set_raster(raster);
}
} }
} }
} }

10
src/image_scaling.cpp
View file

@ -133,6 +133,7 @@ void scale_image_agg(T & target, T const& source, scaling_method_e scaling_metho
// create a scaling matrix // create a scaling matrix
agg::trans_affine img_mtx; agg::trans_affine img_mtx;
img_mtx *= agg::trans_affine_translation(x_off_f, y_off_f);
img_mtx /= agg::trans_affine_scaling(image_ratio_x, image_ratio_y); img_mtx /= agg::trans_affine_scaling(image_ratio_x, image_ratio_y);
// create a linear interpolator for our scaling matrix // create a linear interpolator for our scaling matrix
@ -141,11 +142,10 @@ void scale_image_agg(T & target, T const& source, scaling_method_e scaling_metho
double scaled_width = target.width(); double scaled_width = target.width();
double scaled_height = target.height(); double scaled_height = target.height();
ras.reset(); ras.reset();
ras.move_to_d(x_off_f, y_off_f); ras.move_to_d(0.0, 0.0);
ras.line_to_d(x_off_f + scaled_width, y_off_f); ras.line_to_d(scaled_width, 0.0);
ras.line_to_d(x_off_f + scaled_width, y_off_f + scaled_height); ras.line_to_d(scaled_width, scaled_height);
ras.line_to_d(x_off_f, y_off_f + scaled_height); ras.line_to_d(0.0, scaled_height);
if (scaling_method == SCALING_NEAR) if (scaling_method == SCALING_NEAR)
{ {
using span_gen_type = typename detail::agg_scaling_traits<image_type>::span_image_filter; using span_gen_type = typename detail::agg_scaling_traits<image_type>::span_image_filter;

2
test/data-visual

@ -1 +1 @@
Subproject commit 6393beb4d84448c36b7653abbfabb47e7070bb7c Subproject commit 34648ac6941c09b19919822c867378d1f30915c9