From 170e20e864c2ca8bd1ccc2a94f412bed518d29d4 Mon Sep 17 00:00:00 2001 From: Blake Thompson Date: Mon, 8 May 2017 22:37:35 -0500 Subject: [PATCH 1/7] First attempt at making raster overzooming and scaling work correctly and handle offsets properly. --- include/mapnik/raster.hpp | 13 ++++- .../process_raster_symbolizer.hpp | 48 ++++++++++------- plugins/input/raster/raster_featureset.cpp | 51 ++++++++++--------- src/image_scaling.cpp | 10 ++-- test/data-visual | 2 +- 5 files changed, 76 insertions(+), 48 deletions(-) diff --git a/include/mapnik/raster.hpp b/include/mapnik/raster.hpp index 02022d359..215611946 100644 --- a/include/mapnik/raster.hpp +++ b/include/mapnik/raster.hpp @@ -40,15 +40,27 @@ class raster : private util::noncopyable { public: box2d ext_; + box2d query_ext_; image_any data_; double filter_factor_; boost::optional nodata_; + + template + raster(box2d const& ext, + box2d const& query_ext, + ImageData && data, + double filter_factor) + : ext_(ext), + query_ext_(query_ext), + data_(std::move(data)), + filter_factor_(filter_factor) {} template raster(box2d const& ext, ImageData && data, double filter_factor) : ext_(ext), + query_ext_(ext), data_(std::move(data)), filter_factor_(filter_factor) {} @@ -71,7 +83,6 @@ public: { filter_factor_ = factor; } - }; } diff --git a/include/mapnik/renderer_common/process_raster_symbolizer.hpp b/include/mapnik/renderer_common/process_raster_symbolizer.hpp index 7ca90a0af..9c2853fdf 100644 --- a/include/mapnik/renderer_common/process_raster_symbolizer.hpp +++ b/include/mapnik/renderer_common/process_raster_symbolizer.hpp @@ -53,6 +53,7 @@ struct image_dispatcher image_dispatcher(int start_x, int start_y, int width, int height, double scale_x, double scale_y, + double offset_x, double offset_y, scaling_method_e method, double filter_factor, double opacity, composite_mode_e comp_op, raster_symbolizer const& sym, feature_impl const& feature, @@ -63,15 +64,17 @@ struct image_dispatcher height_(height), scale_x_(scale_x), scale_y_(scale_y), + offset_x_(offset_x), + offset_y_(offset_y), method_(method), filter_factor_(filter_factor), - opacity_(opacity), - comp_op_(comp_op), - sym_(sym), - feature_(feature), - composite_(composite), - nodata_(nodata), - need_scaling_(need_scaling) {} + opacity_(opacity), + comp_op_(comp_op), + sym_(sym), + feature_(feature), + composite_(composite), + nodata_(nodata), + need_scaling_(need_scaling) {} void operator() (image_null const&) const {} //no-op void operator() (image_rgba8 const& data_in) const @@ -79,7 +82,7 @@ struct image_dispatcher if (need_scaling_) { 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_); } else @@ -97,7 +100,7 @@ struct image_dispatcher if (need_scaling_) { 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_); } else @@ -114,6 +117,8 @@ private: int height_; double scale_x_; double scale_y_; + double offset_x_; + double offset_y_; scaling_method_e method_; double filter_factor_; double opacity_; @@ -210,12 +215,17 @@ void render_raster_symbolizer(raster_symbolizer const& sym, if (source) { box2d target_ext = box2d(source->ext_); - prj_trans.backward(target_ext, PROJ_ENVELOPE_POINTS); + box2d target_query_ext = box2d(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 ext = common.t_.forward(target_ext); - int start_x = static_cast(std::floor(ext.minx()+.5)); - int start_y = static_cast(std::floor(ext.miny()+.5)); - int end_x = static_cast(std::floor(ext.maxx()+.5)); - int end_y = static_cast(std::floor(ext.maxy()+.5)); + box2d query_ext = common.t_.forward(target_query_ext); + int start_x = static_cast(std::floor(query_ext.minx()+.5)); + int start_y = static_cast(std::floor(query_ext.miny()+.5)); + int end_x = static_cast(std::floor(query_ext.maxx()+.5)); + int end_y = static_cast(std::floor(query_ext.maxy()+.5)); int raster_width = end_x - start_x; int raster_height = end_y - start_y; if (raster_width > 0 && raster_height > 0) @@ -236,17 +246,20 @@ void render_raster_symbolizer(raster_symbolizer const& sym, if (!prj_trans.equal()) { - double offset_x = ext.minx() - start_x; - double offset_y = ext.miny() - start_y; + // This path does not currently work and is still being figured out. + double offset_x = query_ext.minx() - start_x; + double offset_y = query_ext.miny() - start_y; unsigned mesh_size = static_cast(get(sym,keys::mesh_size,feature, common.vars_, 16)); detail::image_warp_dispatcher 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(), opacity, comp_op, sym, feature, composite, source->nodata()); util::apply_visitor(dispatcher, source->data_); } 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_y = ext.height() / source->data_.height(); double eps = 1e-5; @@ -256,6 +269,7 @@ void render_raster_symbolizer(raster_symbolizer const& sym, (std::abs(start_y) > eps); detail::image_dispatcher dispatcher(start_x, start_y, raster_width, raster_height, image_ratio_x, image_ratio_y, + offset_x, offset_y, scaling_method, source->get_filter_factor(), opacity, comp_op, sym, feature, composite, source->nodata(), scale); util::apply_visitor(dispatcher, source->data_); diff --git a/plugins/input/raster/raster_featureset.cpp b/plugins/input/raster/raster_featureset.cpp index 26cac964f..3228ce529 100644 --- a/plugins/input/raster/raster_featureset.cpp +++ b/plugins/input/raster/raster_featureset.cpp @@ -89,33 +89,36 @@ feature_ptr raster_featureset::next() box2d intersect = bbox_.intersect(curIter_->envelope()); box2d ext = t.forward(intersect); box2d rem = policy_.transform(ext); - if (ext.width() > 0.5 && ext.height() > 0.5 ) - { - // select minimum raster containing whole ext - int x_off = static_cast(std::floor(ext.minx())); - int y_off = static_cast(std::floor(ext.miny())); - int end_x = static_cast(std::ceil(ext.maxx())); - int end_y = static_cast(std::ceil(ext.maxy())); + // select minimum raster containing whole ext + int x_off = static_cast(std::floor(ext.minx())); + int y_off = static_cast(std::floor(ext.miny())); + int end_x = static_cast(std::ceil(ext.maxx())); + int end_y = static_cast(std::ceil(ext.maxy())); - // clip to available data - if (x_off < 0) x_off = 0; - 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; + // clip to available data + if (x_off < 0) x_off = 0; + 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 height = end_y - y_off; - - // calculate actual box2d of returned raster - box2d feature_raster_extent(rem.minx() + x_off, - rem.miny() + y_off, - rem.maxx() + x_off + width, - rem.maxy() + y_off + height); - intersect = t.backward(feature_raster_extent); - mapnik::image_any data = reader->read(x_off, y_off, width, height); - mapnik::raster_ptr raster = std::make_shared(intersect, std::move(data), 1.0); - feature->set_raster(raster); + int width = end_x - x_off; + int height = end_y - y_off; + if (width < 1) { + width = 1; } + if (height < 1) { + height = 1; + } + + // calculate actual box2d of returned raster + box2d feature_raster_extent(rem.minx() + x_off, + rem.miny() + y_off, + rem.maxx() + x_off + width, + rem.maxy() + y_off + height); + feature_raster_extent = t.backward(feature_raster_extent); + mapnik::image_any data = reader->read(x_off, y_off, width, height); + mapnik::raster_ptr raster = std::make_shared(feature_raster_extent, intersect, std::move(data), 1.0); + feature->set_raster(raster); } } } diff --git a/src/image_scaling.cpp b/src/image_scaling.cpp index 413948c64..4c6bb49c1 100644 --- a/src/image_scaling.cpp +++ b/src/image_scaling.cpp @@ -133,6 +133,7 @@ void scale_image_agg(T & target, T const& source, scaling_method_e scaling_metho // create a scaling matrix 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); // 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_height = target.height(); ras.reset(); - ras.move_to_d(x_off_f, y_off_f); - ras.line_to_d(x_off_f + scaled_width, y_off_f); - ras.line_to_d(x_off_f + scaled_width, y_off_f + scaled_height); - ras.line_to_d(x_off_f, y_off_f + scaled_height); - + ras.move_to_d(0.0, 0.0); + ras.line_to_d(scaled_width, 0.0); + ras.line_to_d(scaled_width, scaled_height); + ras.line_to_d(0.0, scaled_height); if (scaling_method == SCALING_NEAR) { using span_gen_type = typename detail::agg_scaling_traits::span_image_filter; diff --git a/test/data-visual b/test/data-visual index 17a3e7122..d96d10677 160000 --- a/test/data-visual +++ b/test/data-visual @@ -1 +1 @@ -Subproject commit 17a3e712266a6fac4d89c8473fc0429f6c54fae3 +Subproject commit d96d1067796f1fe99da27941b216cb6991661a6a From 85d7f286103854864ad9d913b02499f84fae588b Mon Sep 17 00:00:00 2001 From: artemp Date: Tue, 9 May 2017 09:54:51 +0200 Subject: [PATCH 2/7] remove QMAKE_MAC_SDK requirement - just use what is available --- demo/viewer/viewer.pro | 1 - 1 file changed, 1 deletion(-) diff --git a/demo/viewer/viewer.pro b/demo/viewer/viewer.pro index 8a9747752..aad00b424 100644 --- a/demo/viewer/viewer.pro +++ b/demo/viewer/viewer.pro @@ -1,7 +1,6 @@ ###################################################################### # Mapnik viewer - Copyright (C) 2007 Artem Pavlenko ###################################################################### -QMAKE_MAC_SDK = macosx10.11 TEMPLATE = app QT += core gui widgets QMAKE_CXX = $$system(mapnik-config --cxx) From adacb16a2c4607454ffd3582e34b634a53f496be Mon Sep 17 00:00:00 2001 From: Blake Thompson Date: Tue, 9 May 2017 12:15:03 -0500 Subject: [PATCH 3/7] Removed 0.5 limit on width and height in gdal plugin, added new visual tests --- plugins/input/gdal/gdal_featureset.cpp | 10 ---------- test/data-visual | 2 +- 2 files changed, 1 insertion(+), 11 deletions(-) diff --git a/plugins/input/gdal/gdal_featureset.cpp b/plugins/input/gdal/gdal_featureset.cpp index 1377322e3..32ebc4863 100644 --- a/plugins/input/gdal/gdal_featureset.cpp +++ b/plugins/input/gdal/gdal_featureset.cpp @@ -187,16 +187,6 @@ 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 feature_raster_extent(x_off, y_off, x_off + width, y_off + height); intersect = t.backward(feature_raster_extent); diff --git a/test/data-visual b/test/data-visual index d96d10677..457eb2194 160000 --- a/test/data-visual +++ b/test/data-visual @@ -1 +1 @@ -Subproject commit d96d1067796f1fe99da27941b216cb6991661a6a +Subproject commit 457eb219475d5a50e742d97ccbbe1015ff4f6246 From 318be793f92444ea1f8c07002954af3302aa22c1 Mon Sep 17 00:00:00 2001 From: Blake Thompson Date: Wed, 10 May 2017 09:04:29 -0500 Subject: [PATCH 4/7] Removed filter factor --- .../process_raster_symbolizer.hpp | 3 +- plugins/input/gdal/gdal_featureset.cpp | 710 +++++++++--------- plugins/input/raster/raster_featureset.cpp | 10 +- 3 files changed, 351 insertions(+), 372 deletions(-) diff --git a/include/mapnik/renderer_common/process_raster_symbolizer.hpp b/include/mapnik/renderer_common/process_raster_symbolizer.hpp index 9c2853fdf..54726eaf5 100644 --- a/include/mapnik/renderer_common/process_raster_symbolizer.hpp +++ b/include/mapnik/renderer_common/process_raster_symbolizer.hpp @@ -216,7 +216,8 @@ void render_raster_symbolizer(raster_symbolizer const& sym, { box2d target_ext = box2d(source->ext_); box2d target_query_ext = box2d(source->query_ext_); - if (!prj_trans.equal()) { + if (!prj_trans.equal()) + { prj_trans.backward(target_ext, PROJ_ENVELOPE_POINTS); prj_trans.backward(target_query_ext, PROJ_ENVELOPE_POINTS); } diff --git a/plugins/input/gdal/gdal_featureset.cpp b/plugins/input/gdal/gdal_featureset.cpp index 32ebc4863..f71535fd6 100644 --- a/plugins/input/gdal/gdal_featureset.cpp +++ b/plugins/input/gdal/gdal_featureset.cpp @@ -189,7 +189,7 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q) //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); + 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; @@ -198,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::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(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::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(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::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::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(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::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(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::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 << "," <c2 << "," << ce->c3; + } } - mapnik::raster_ptr raster = std::make_shared(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::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(&imageData[i]) = 0; + } + else + { + *reinterpret_cast(&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(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::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(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::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(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::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 << "," <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(&imageData[i]) = 0; } else { - MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming gray band)"; - grey = band; + *reinterpret_cast(&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(&imageData[i]) = 0; + row[x] = (ce->c4 << 24)| (ce->c3 << 16) | (ce->c2 << 8) | (ce->c1) ; } else { - *reinterpret_cast(&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(&imageData[i]) = 0; - } - else - { - *reinterpret_cast(&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()); @@ -573,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(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(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(); } diff --git a/plugins/input/raster/raster_featureset.cpp b/plugins/input/raster/raster_featureset.cpp index 3228ce529..edea94322 100644 --- a/plugins/input/raster/raster_featureset.cpp +++ b/plugins/input/raster/raster_featureset.cpp @@ -96,6 +96,8 @@ feature_ptr raster_featureset::next() int end_y = static_cast(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::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 feature_raster_extent(rem.minx() + x_off, From bc9dcdc5843ecdcfaa52229f74cdaa4e9d4568ca Mon Sep 17 00:00:00 2001 From: Blake Thompson Date: Wed, 10 May 2017 13:31:11 -0500 Subject: [PATCH 5/7] Updated visual test data --- test/data-visual | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/data-visual b/test/data-visual index 457eb2194..7023ba0a7 160000 --- a/test/data-visual +++ b/test/data-visual @@ -1 +1 @@ -Subproject commit 457eb219475d5a50e742d97ccbbe1015ff4f6246 +Subproject commit 7023ba0a782ed64ebe254e653e0f7a9a7ee1d92b From baf6c9f3210129da190bc502792d80365871033c Mon Sep 17 00:00:00 2001 From: Blake Thompson Date: Thu, 11 May 2017 09:57:18 -0500 Subject: [PATCH 6/7] Updated tests, changed boostrap to use latest released mason --- bootstrap.sh | 2 +- test/data-visual | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/bootstrap.sh b/bootstrap.sh index b2b171f49..f97d598f7 100755 --- a/bootstrap.sh +++ b/bootstrap.sh @@ -11,7 +11,7 @@ todo - shrink icu data ' -MASON_VERSION="3c6df04" +MASON_VERSION="v0.10.0" function setup_mason() { if [[ ! -d ./.mason ]]; then diff --git a/test/data-visual b/test/data-visual index 7023ba0a7..d9ba85931 160000 --- a/test/data-visual +++ b/test/data-visual @@ -1 +1 @@ -Subproject commit 7023ba0a782ed64ebe254e653e0f7a9a7ee1d92b +Subproject commit d9ba85931081cc0eec131c6ae654afd54c5798e2 From 7222129e6cf16a7524673037405698796d09b50f Mon Sep 17 00:00:00 2001 From: Blake Thompson Date: Fri, 12 May 2017 10:55:48 -0500 Subject: [PATCH 7/7] Update visual test data to 3.0.x branch --- test/data-visual | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/data-visual b/test/data-visual index d9ba85931..34648ac69 160000 --- a/test/data-visual +++ b/test/data-visual @@ -1 +1 @@ -Subproject commit d9ba85931081cc0eec131c6ae654afd54c5798e2 +Subproject commit 34648ac6941c09b19919822c867378d1f30915c9