support raster reprojection - adapted from work at https://github.com/albertov/mapnik2/ - nice work alberto and meteogrid - closes #663
This commit is contained in:
Dane Springmeyer
parent
67df69834f
commit
9bda9c7019
3 changed files with 122 additions and 28 deletions
|
|
@ -24,29 +24,95 @@
|
|||
// mapnik
|
||||
#include <mapnik/agg_renderer.hpp>
|
||||
#include <mapnik/agg_rasterizer.hpp>
|
||||
#include <mapnik/config.hpp>
|
||||
#include <mapnik/image_data.hpp>
|
||||
#include <mapnik/image_util.hpp>
|
||||
#include <mapnik/raster.hpp>
|
||||
#include <mapnik/box2d.hpp>
|
||||
|
||||
// stl
|
||||
#include <cmath>
|
||||
|
||||
namespace mapnik {
|
||||
|
||||
static inline void resample_raster(raster &target, raster const& source,
|
||||
proj_transform const& prj_trans,
|
||||
double offset_x, double offset_y,
|
||||
double filter_factor,
|
||||
std::string scaling_method_name)
|
||||
{
|
||||
if (prj_trans.equal()) {
|
||||
|
||||
if (scaling_method_name == "bilinear8"){
|
||||
scale_image_bilinear8<image_data_32>(target.data_,source.data_,
|
||||
offset_x, offset_y);
|
||||
} else {
|
||||
double scale_factor = target.data_.width() / source.data_.width();
|
||||
scaling_method_e scaling_method = get_scaling_method_by_name(scaling_method_name);
|
||||
scale_image_agg<image_data_32>(target.data_,source.data_, (scaling_method_e)scaling_method, scale_factor, offset_x, offset_y, filter_factor);
|
||||
}
|
||||
} else {
|
||||
CoordTransform ts(source.data_.width(), source.data_.height(),
|
||||
source.ext_);
|
||||
CoordTransform tt(target.data_.width(), target.data_.height(),
|
||||
target.ext_, offset_x, offset_y);
|
||||
unsigned i, j, mesh_size=16;
|
||||
unsigned mesh_nx = ceil(target.data_.width()/double(mesh_size)+1);
|
||||
unsigned mesh_ny = ceil(target.data_.height()/double(mesh_size)+1);
|
||||
ImageData<double> xs(mesh_nx, mesh_ny);
|
||||
ImageData<double> ys(mesh_nx, mesh_ny);
|
||||
|
||||
// Precalculate reprojected mesh
|
||||
for(j=0; j<mesh_ny; j++) {
|
||||
for (i=0; i<mesh_nx; i++) {
|
||||
xs(i,j) = i*mesh_size+.5;
|
||||
ys(i,j) = j*mesh_size+.5;
|
||||
tt.backward(&xs(i,j), &ys(i,j));
|
||||
}
|
||||
}
|
||||
prj_trans.forward(xs.getData(), ys.getData(), NULL, mesh_nx*mesh_ny);
|
||||
|
||||
// Interpolate. TODO: Use AGG for this
|
||||
for(j=0; j<mesh_ny-1; j++) {
|
||||
for (i=0; i<mesh_nx-1; i++) {
|
||||
box2d<double> win(xs(i,j), ys(i,j), xs(i+1,j+1), ys(i+1,j+1));
|
||||
win = ts.forward(win);
|
||||
CoordTransform tw(mesh_size, mesh_size, win);
|
||||
unsigned x, y;
|
||||
for (y=0; y<mesh_size; y++) {
|
||||
for (x=0; x<mesh_size; x++) {
|
||||
double x2=x, y2=y;
|
||||
tw.backward(&x2,&y2);
|
||||
if (x2>=0 && x2<source.data_.width() &&
|
||||
y2>=0 && y2<source.data_.height())
|
||||
{
|
||||
target.data_(i*mesh_size+x,(j+1)*mesh_size-y) =\
|
||||
source.data_((unsigned)x2,(unsigned)y2);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void agg_renderer<T>::process(raster_symbolizer const& sym,
|
||||
Feature const& feature,
|
||||
proj_transform const& prj_trans)
|
||||
{
|
||||
raster_ptr const& raster=feature.get_raster();
|
||||
if (raster)
|
||||
raster_ptr const& source=feature.get_raster();
|
||||
if (source)
|
||||
{
|
||||
// If there's a colorizer defined, use it to color the raster in-place
|
||||
raster_colorizer_ptr colorizer = sym.get_colorizer();
|
||||
if (colorizer)
|
||||
colorizer->colorize(raster,feature.props());
|
||||
|
||||
box2d<double> ext=t_.forward(raster->ext_);
|
||||
|
||||
colorizer->colorize(source,feature.props());
|
||||
|
||||
box2d<double> target_ext = box2d<double>(source->ext_);
|
||||
prj_trans.backward(target_ext, 20);
|
||||
|
||||
box2d<double> ext=t_.forward(target_ext);
|
||||
int start_x = (int)ext.minx();
|
||||
int start_y = (int)ext.miny();
|
||||
int end_x = (int)ceil(ext.maxx());
|
||||
|
|
@ -56,45 +122,42 @@ void agg_renderer<T>::process(raster_symbolizer const& sym,
|
|||
double err_offs_x = ext.minx() - start_x;
|
||||
double err_offs_y = ext.miny() - start_y;
|
||||
|
||||
if ( raster_width > 0 && raster_height > 0)
|
||||
if (raster_width > 0 && raster_height > 0)
|
||||
{
|
||||
double scale_factor = ext.width() / raster->data_.width();
|
||||
image_data_32 target(raster_width,raster_height);
|
||||
|
||||
if (sym.get_scaling() == "bilinear8"){
|
||||
scale_image_bilinear8<image_data_32>(target,raster->data_, err_offs_x, err_offs_y);
|
||||
} else {
|
||||
scaling_method_e scaling_method = get_scaling_method_by_name(sym.get_scaling());
|
||||
scale_image_agg<image_data_32>(target,raster->data_, (scaling_method_e)scaling_method, scale_factor, err_offs_x, err_offs_y, sym.calculate_filter_factor());
|
||||
}
|
||||
image_data_32 target_data(raster_width,raster_height);
|
||||
raster target(target_ext, target_data);
|
||||
|
||||
resample_raster(target, *source, prj_trans, err_offs_x, err_offs_y,
|
||||
sym.calculate_filter_factor(),
|
||||
sym.get_scaling());
|
||||
|
||||
if (sym.get_mode() == "normal"){
|
||||
if (sym.get_opacity() == 1.0) {
|
||||
pixmap_.set_rectangle_alpha(start_x,start_y,target);
|
||||
pixmap_.set_rectangle_alpha(start_x,start_y,target.data_);
|
||||
} else {
|
||||
pixmap_.set_rectangle_alpha2(target,start_x,start_y, sym.get_opacity());
|
||||
pixmap_.set_rectangle_alpha2(target.data_,start_x,start_y, sym.get_opacity());
|
||||
}
|
||||
} else if (sym.get_mode() == "grain_merge"){
|
||||
pixmap_.template merge_rectangle<MergeGrain> (target,start_x,start_y, sym.get_opacity());
|
||||
pixmap_.template merge_rectangle<MergeGrain> (target.data_,start_x,start_y, sym.get_opacity());
|
||||
} else if (sym.get_mode() == "grain_merge2"){
|
||||
pixmap_.template merge_rectangle<MergeGrain2> (target,start_x,start_y, sym.get_opacity());
|
||||
pixmap_.template merge_rectangle<MergeGrain2> (target.data_,start_x,start_y, sym.get_opacity());
|
||||
} else if (sym.get_mode() == "multiply"){
|
||||
pixmap_.template merge_rectangle<Multiply> (target,start_x,start_y, sym.get_opacity());
|
||||
pixmap_.template merge_rectangle<Multiply> (target.data_,start_x,start_y, sym.get_opacity());
|
||||
} else if (sym.get_mode() == "multiply2"){
|
||||
pixmap_.template merge_rectangle<Multiply2> (target,start_x,start_y, sym.get_opacity());
|
||||
pixmap_.template merge_rectangle<Multiply2> (target.data_,start_x,start_y, sym.get_opacity());
|
||||
} else if (sym.get_mode() == "divide"){
|
||||
pixmap_.template merge_rectangle<Divide> (target,start_x,start_y, sym.get_opacity());
|
||||
pixmap_.template merge_rectangle<Divide> (target.data_,start_x,start_y, sym.get_opacity());
|
||||
} else if (sym.get_mode() == "divide2"){
|
||||
pixmap_.template merge_rectangle<Divide2> (target,start_x,start_y, sym.get_opacity());
|
||||
pixmap_.template merge_rectangle<Divide2> (target.data_,start_x,start_y, sym.get_opacity());
|
||||
} else if (sym.get_mode() == "screen"){
|
||||
pixmap_.template merge_rectangle<Screen> (target,start_x,start_y, sym.get_opacity());
|
||||
pixmap_.template merge_rectangle<Screen> (target.data_,start_x,start_y, sym.get_opacity());
|
||||
} else if (sym.get_mode() == "hard_light"){
|
||||
pixmap_.template merge_rectangle<HardLight> (target,start_x,start_y, sym.get_opacity());
|
||||
pixmap_.template merge_rectangle<HardLight> (target.data_,start_x,start_y, sym.get_opacity());
|
||||
} else {
|
||||
if (sym.get_opacity() == 1.0){
|
||||
pixmap_.set_rectangle(start_x,start_y,target);
|
||||
pixmap_.set_rectangle(start_x,start_y,target.data_);
|
||||
} else {
|
||||
pixmap_.set_rectangle_alpha2(target,start_x,start_y, sym.get_opacity());
|
||||
pixmap_.set_rectangle_alpha2(target.data_,start_x,start_y, sym.get_opacity());
|
||||
}
|
||||
}
|
||||
// TODO: other modes? (add,diff,sub,...)
|
||||
|
|
|
|||
|
|
@ -1373,6 +1373,7 @@ void cairo_renderer_base::process(polygon_pattern_symbolizer const& sym,
|
|||
}
|
||||
}
|
||||
|
||||
//FIXME: Port reprojection code from agg/process_raster_symbolizer.cpp
|
||||
void cairo_renderer_base::process(raster_symbolizer const& sym,
|
||||
Feature const& feature,
|
||||
proj_transform const& /*prj_trans*/)
|
||||
|
|
|
|||
|
|
@ -131,6 +131,36 @@ def test_raster_with_alpha_blends_correctly_with_background():
|
|||
imdata = mim.tostring()
|
||||
# All white is expected
|
||||
assert contains_word('\xff\xff\xff\xff', imdata)
|
||||
|
||||
def test_raster_warping():
|
||||
lyrSrs = "+init=epsg:32630"
|
||||
mapSrs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
|
||||
lyr = mapnik2.Layer('dataraster', lyrSrs)
|
||||
lyr.datasource = mapnik2.Gdal(
|
||||
file = '../data/raster/dataraster.tif',
|
||||
band = 1,
|
||||
)
|
||||
sym = mapnik2.RasterSymbolizer()
|
||||
sym.colorizer = mapnik2.RasterColorizer(mapnik2.COLORIZER_DISCRETE, mapnik2.Color(255,255,0))
|
||||
rule = mapnik2.Rule()
|
||||
rule.symbols.append(sym)
|
||||
style = mapnik2.Style()
|
||||
style.rules.append(rule)
|
||||
_map = mapnik2.Map(256,256, mapSrs)
|
||||
_map.append_style('foo', style)
|
||||
lyr.styles.append('foo')
|
||||
_map.layers.append(lyr)
|
||||
prj_trans = mapnik2.ProjTransform(mapnik2.Projection(mapSrs),
|
||||
mapnik2.Projection(lyrSrs))
|
||||
_map.zoom_to_box(prj_trans.backward(lyr.envelope()))
|
||||
|
||||
im = mapnik2.Image(_map.width,_map.height)
|
||||
mapnik2.render(_map, im)
|
||||
#print mapnik2.save_map_to_string(_map)
|
||||
# save a png somewhere so we can see it
|
||||
save_data('test_raster_warping.png', im.tostring('png'))
|
||||
imdata = im.tostring()
|
||||
assert contains_word('\xff\xff\x00\xff', imdata)
|
||||
|
||||
if __name__ == "__main__":
|
||||
setup()
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue