Implement input plugin for PostGIS Rasters

Closes #1660.
Includes python bindings and testcases.
This commit is contained in:
Sandro Santilli 2014-06-26 12:38:17 +02:00
parent 916a30124b
commit fa4d1e4ada
12 changed files with 3278 additions and 1 deletions

View file

@ -100,6 +100,7 @@ pretty_dep_names = {
PLUGINS = { # plugins with external dependencies
# configured by calling project, hence 'path':None
'postgis': {'default':True,'path':None,'inc':'libpq-fe.h','lib':'pq','lang':'C'},
'pgraster': {'default':True,'path':None,'inc':'libpq-fe.h','lib':'pq','lang':'C'},
'gdal': {'default':True,'path':None,'inc':'gdal_priv.h','lib':'gdal','lang':'C++'},
'ogr': {'default':True,'path':None,'inc':'ogrsf_frmts.h','lib':'gdal','lang':'C++'},
# configured with custom paths, hence 'path': PREFIX/INCLUDES/LIBS
@ -1381,7 +1382,7 @@ if not preconfigured:
env['LIBS'].remove(libname)
else:
details['lib'] = libname
elif plugin == 'postgis':
elif plugin == 'postgis' or plugin == 'pgraster':
conf.parse_pg_config('PG_CONFIG')
elif plugin == 'ogr':
if conf.ogr_enabled():

View file

@ -439,6 +439,52 @@ def PostGIS(**keywords):
keywords['type'] = 'postgis'
return CreateDatasource(keywords)
def PgRaster(**keywords):
"""Create a PgRaster Datasource.
Required keyword arguments:
dbname -- database name to connect to
table -- table name or subselect query
*Note: if using subselects for the 'table' value consider also
passing the 'raster_field' and 'srid' and 'extent_from_subquery'
options and/or specifying the 'raster_table' option.
Optional db connection keyword arguments:
user -- database user to connect as (default: see postgres docs)
password -- password for database user (default: see postgres docs)
host -- portgres hostname (default: see postgres docs)
port -- postgres port (default: see postgres docs)
initial_size -- integer size of connection pool (default: 1)
max_size -- integer max of connection pool (default: 10)
persist_connection -- keep connection open (default: True)
Optional table-level keyword arguments:
extent -- manually specified data extent (comma delimited string, default: None)
estimate_extent -- boolean, direct PostGIS to use the faster, less accurate `estimate_extent` over `extent` (default: False)
extent_from_subquery -- boolean, direct Mapnik to query Postgis for the extent of the raw 'table' value (default: uses 'geometry_table')
raster_table -- specify geometry table to use to look up metadata (default: automatically parsed from 'table' value)
raster_field -- specify geometry field to use (default: first entry in raster_columns)
srid -- specify srid to use (default: auto-detected from geometry_field)
row_limit -- integer limit of rows to return (default: 0)
cursor_size -- integer size of binary cursor to use (default: 0, no binary cursor is used)
use_overviews -- boolean, use overviews when available (default: false)
prescale_rasters -- boolean, scale rasters on the db side (default: false)
clip_rasters -- boolean, clip rasters on the db side (default: false)
band -- integer, if non-zero interprets the given band (1-based offset) as a data raster (default: 0)
>>> from mapnik import PgRaster, Layer
>>> params = dict(dbname='mapnik',table='osm',user='postgres',password='gis')
>>> params['estimate_extent'] = False
>>> params['extent'] = '-20037508,-19929239,20037508,19929239'
>>> pgraster = PgRaster(**params)
>>> lyr = Layer('PgRaster Layer')
>>> lyr.datasource = pgraster
"""
keywords['type'] = 'pgraster'
return CreateDatasource(keywords)
def Raster(**keywords):
"""Create a Raster (Tiff) Datasource.

View file

@ -0,0 +1,14 @@
This plugin shares directives with the "postgis" one,
with the following differences:
- "raster_field" replaces "geometry_field"
- "raster_table" replaces "geometry_table"
- "prescale_rasters" replaces "simplify_geometries"
- "use_overviews" introduced
- "clip_rasters" boolean introduced, defaults to false
- "band" introduced, with same semantic of the GDAL driver

View file

@ -0,0 +1,23 @@
- Automated tests
- Test subqueries [x]
- Test Data input [x]
- Test RGB input [x]
- Test RGBA input [x]
- Test Grayscale input [x]
- Support for all band types:
- PT_1BB data[x] rgb[ ] grayscale[x]
- PT_2BUI data[x] rgb[ ] grayscale[x]
- PT_4BUI data[x] rgb[ ] grayscale[x]
- PT_8BSI data[x] rgb[x] grayscale[x]
- PT_8BUI data[x] rgb[x] grayscale[x]
- PT_16BSI data[x] rgb[ ] grayscale[x]
- PT_16BUI data[x] rgb[ ] grayscale[x]
- PT_32BSI data[x] rgb[ ] grayscale[x]
- PT_32BUI data[x] rgb[ ] grayscale[x]
- PT_32BF data[x] rgb[ ] grayscale[ ]
- PT_64BF data[x] rgb[ ] grayscale[ ]
- Have pgraster and postgis plugins share the same connection pool
- Make overviews enabled by default if table is not a subquery ?
- Make clipping enabled automatically when needed ?
- Allow more flexible band layout specification, see
http://github.com/mapnik/mapnik/wiki/RFC:-Raster-color-interpretation

View file

@ -0,0 +1,84 @@
#
# This file is part of Mapnik (c++ mapping toolkit)
#
# Copyright (C) 2013 Artem Pavlenko
#
# Mapnik is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with this library; if not, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
#
#
Import ('plugin_base')
Import ('env')
from copy import copy
PLUGIN_NAME = 'pgraster'
plugin_env = plugin_base.Clone()
plugin_sources = Split(
"""
%(PLUGIN_NAME)s_datasource.cpp
%(PLUGIN_NAME)s_featureset.cpp
%(PLUGIN_NAME)s_wkb_reader.cpp
""" % locals()
)
cxxflags = []
plugin_env['LIBS'] = []
# TODO: use variable for top-level dir
plugin_env['CXXFLAGS'].append('-Iplugins/input/postgis')
if env['RUNTIME_LINK'] == 'static':
# pkg-config is more reliable than pg_config across platforms
cmd = 'pkg-config libpq --libs --static'
try:
plugin_env.ParseConfig(cmd)
except OSError, e:
plugin_env.Append(LIBS='pq')
else:
plugin_env.Append(LIBS='pq')
# Link Library to Dependencies
libraries = copy(plugin_env['LIBS'])
if env['THREADING'] == 'multi':
libraries.append('boost_thread%s' % env['BOOST_APPEND'])
if env['PLUGIN_LINKING'] == 'shared':
libraries.insert(0,env['MAPNIK_NAME'])
libraries.append(env['ICU_LIB_NAME'])
libraries.append('boost_system%s' % env['BOOST_APPEND'])
TARGET = plugin_env.SharedLibrary('../%s' % PLUGIN_NAME,
SHLIBPREFIX='',
SHLIBSUFFIX='.input',
source=plugin_sources,
LIBS=libraries,
LINKFLAGS=env['CUSTOM_LDFLAGS'])
# if the plugin links to libmapnik ensure it is built first
Depends(TARGET, env.subst('../../../src/%s' % env['MAPNIK_LIB_NAME']))
if 'uninstall' not in COMMAND_LINE_TARGETS:
env.Install(env['MAPNIK_INPUT_PLUGINS_DEST'], TARGET)
env.Alias('install', env['MAPNIK_INPUT_PLUGINS_DEST'])
plugin_obj = {
'LIBS': libraries,
'SOURCES': plugin_sources,
}
Return('plugin_obj')

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,147 @@
/*****************************************************************************
*
* This file is part of Mapnik (c++ mapping toolkit)
*
* Copyright (C) 2014 Artem Pavlenko
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
*****************************************************************************
*
* Initially developed by Sandro Santilli <strk@keybit.net>
*
*****************************************************************************/
#ifndef PGRASTER_DATASOURCE_HPP
#define PGRASTER_DATASOURCE_HPP
// mapnik
#include <mapnik/datasource.hpp>
#include <mapnik/params.hpp>
#include <mapnik/query.hpp>
#include <mapnik/feature.hpp>
#include <mapnik/box2d.hpp>
#include <mapnik/coord.hpp>
#include <mapnik/feature_layer_desc.hpp>
#include <mapnik/unicode.hpp>
#include <mapnik/value_types.hpp>
// boost
#include <boost/optional.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/scoped_ptr.hpp>
// stl
#include <vector>
#include <string>
#include "connection_manager.hpp"
#include "resultset.hpp"
#include "cursorresultset.hpp"
using mapnik::transcoder;
using mapnik::datasource;
using mapnik::feature_style_context_map;
using mapnik::processor_context_ptr;
using mapnik::box2d;
using mapnik::layer_descriptor;
using mapnik::featureset_ptr;
using mapnik::feature_ptr;
using mapnik::query;
using mapnik::parameters;
using mapnik::coord2d;
typedef boost::shared_ptr< ConnectionManager::PoolType> CnxPool_ptr;
struct pgraster_overview
{
std::string schema;
std::string table;
std::string column;
float scale; // max absolute scale between x and y
};
class pgraster_datasource : public datasource
{
public:
pgraster_datasource(const parameters &params);
~pgraster_datasource();
mapnik::datasource::datasource_t type() const;
static const char * name();
processor_context_ptr get_context(feature_style_context_map &) const;
featureset_ptr features_with_context(query const& q, processor_context_ptr ctx) const;
featureset_ptr features(query const& q) const;
featureset_ptr features_at_point(coord2d const& pt, double tol = 0) const;
mapnik::box2d<double> envelope() const;
boost::optional<mapnik::datasource::geometry_t> get_geometry_type() const;
layer_descriptor get_descriptor() const;
private:
std::string sql_bbox(box2d<double> const& env) const;
std::string populate_tokens(std::string const& sql, double scale_denom, box2d<double> const& env, double pixel_width, double pixel_height) const;
std::string populate_tokens(std::string const& sql) const;
boost::shared_ptr<IResultSet> get_resultset(boost::shared_ptr<Connection> &conn, std::string const& sql, CnxPool_ptr const& pool, processor_context_ptr ctx= processor_context_ptr()) const;
static const std::string RASTER_COLUMNS;
static const std::string RASTER_OVERVIEWS;
static const std::string SPATIAL_REF_SYS;
static const double FMAX;
const std::string uri_;
const std::string username_;
const std::string password_;
// table name (schema qualified or not) or subquery
const std::string table_;
// schema name (possibly extracted from table_)
std::string schema_;
// table name (possibly extracted from table_)
std::string raster_table_;
const std::string raster_field_;
std::string key_field_;
mapnik::value_integer cursor_fetch_size_;
mapnik::value_integer row_limit_;
std::string geometryColumn_;
mapnik::datasource::datasource_t type_;
int srid_;
// 1-based index of band to extract from the raster
// 0 means fetch all bands
// any index also forces color interpretation off so that values
// arrives untouched into the resulting mapnik raster, for threatment
// by raster colorizer
int band_;
// Available overviews, ordered by max scale, ascending
std::vector<pgraster_overview> overviews_;
mutable bool extent_initialized_;
mutable mapnik::box2d<double> extent_;
bool prescale_rasters_;
bool use_overviews_;
bool clip_rasters_;
layer_descriptor desc_;
ConnectionCreator<Connection> creator_;
const std::string bbox_token_;
const std::string scale_denom_token_;
const std::string pixel_width_token_;
const std::string pixel_height_token_;
int pool_max_size_;
bool persist_connection_;
bool extent_from_subquery_;
bool estimate_extent_;
int max_async_connections_;
bool asynchronous_request_;
int intersect_min_scale_;
int intersect_max_scale_;
};
#endif // PGRASTER_DATASOURCE_HPP

View file

@ -0,0 +1,360 @@
/*****************************************************************************
*
* This file is part of Mapnik (c++ mapping toolkit)
*
* Copyright (C) 2014 Artem Pavlenko
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
*****************************************************************************
*
* Initially developed by Sandro Santilli <strk@keybit.net>
*
*****************************************************************************/
#include "pgraster_featureset.hpp"
#include "pgraster_wkb_reader.hpp"
#include "resultset.hpp"
#include "cursorresultset.hpp"
// mapnik
#include <mapnik/global.hpp>
#include <mapnik/debug.hpp>
#include <mapnik/unicode.hpp>
#include <mapnik/value_types.hpp>
#include <mapnik/feature_factory.hpp>
#include <mapnik/ctrans.hpp>
#include <mapnik/raster.hpp>
#include <mapnik/image_data.hpp>
#include <mapnik/util/conversions.hpp>
#include <mapnik/util/trim.hpp>
#include <mapnik/global.hpp> // for int2net
#include <boost/scoped_array.hpp>
// boost
#include <boost/cstdint.hpp> // for boost::int16_t
// stl
#include <sstream>
#include <string>
using mapnik::geometry_type;
using mapnik::byte;
using mapnik::feature_factory;
using mapnik::context_ptr;
pgraster_featureset::pgraster_featureset(boost::shared_ptr<IResultSet> const& rs,
context_ptr const& ctx,
std::string const& encoding,
bool key_field, int bandno)
: rs_(rs),
ctx_(ctx),
tr_(new transcoder(encoding)),
feature_id_(1),
key_field_(key_field),
band_(bandno)
{
}
std::string numeric2string(const char* buf);
feature_ptr pgraster_featureset::next()
{
while (rs_->next())
{
// new feature
unsigned pos = 1;
feature_ptr feature;
if (key_field_)
{
std::string name = rs_->getFieldName(pos);
// null feature id is not acceptable
if (rs_->isNull(pos))
{
MAPNIK_LOG_WARN(pgraster) << "pgraster_featureset: null value encountered for key_field: " << name;
continue;
}
// create feature with user driven id from attribute
int oid = rs_->getTypeOID(pos);
const char* buf = rs_->getValue(pos);
// validation happens of this type at initialization
mapnik::value_integer val;
if (oid == 20)
{
val = int8net(buf);
}
else if (oid == 21)
{
val = int2net(buf);
}
else
{
val = int4net(buf);
}
MAPNIK_LOG_WARN(pgraster) << "pgraster_featureset: feature key: " << val;
feature = feature_factory::create(ctx_, val);
// TODO - extend feature class to know
// that its id is also an attribute to avoid
// this duplication
feature->put<mapnik::value_integer>(name,val);
++pos;
}
else
{
// fallback to auto-incrementing id
MAPNIK_LOG_WARN(pgraster) << "pgraster_featureset: feature id: " << feature_id_;
feature = feature_factory::create(ctx_, feature_id_);
++feature_id_;
}
// null geometry is not acceptable
if (rs_->isNull(0))
{
MAPNIK_LOG_WARN(pgraster) << "pgraster_featureset: null value encountered for raster";
continue;
}
// parse geometry
int size = rs_->getFieldLength(0);
const uint8_t *data = (const uint8_t*)rs_->getValue(0);
mapnik::raster_ptr raster = pgraster_wkb_reader::read(data, size, band_);
if (!raster)
{
MAPNIK_LOG_WARN(pgraster) << "pgraster_featureset: could not parse raster wkb";
// TODO: throw an exception ?
continue;
}
MAPNIK_LOG_WARN(pgraster) << "pgraster_featureset: raster of " << raster->data_.width() << "x" << raster->data_.height() << " pixels covering extent " << raster->ext_;
feature->set_raster(raster);
unsigned num_attrs = ctx_->size() + 1;
for (; pos < num_attrs; ++pos)
{
std::string name = rs_->getFieldName(pos);
// NOTE: we intentionally do not store null here
// since it is equivalent to the attribute not existing
if (!rs_->isNull(pos))
{
const char* buf = rs_->getValue(pos);
const int oid = rs_->getTypeOID(pos);
switch (oid)
{
case 16: //bool
{
feature->put(name, (buf[0] != 0));
break;
}
case 23: //int4
{
feature->put<mapnik::value_integer>(name, int4net(buf));
break;
}
case 21: //int2
{
feature->put<mapnik::value_integer>(name, int2net(buf));
break;
}
case 20: //int8/BigInt
{
feature->put<mapnik::value_integer>(name, int8net(buf));
break;
}
case 700: //float4
{
float val;
float4net(val, buf);
feature->put(name, static_cast<double>(val));
break;
}
case 701: //float8
{
double val;
float8net(val, buf);
feature->put(name, val);
break;
}
case 25: //text
case 1043: //varchar
case 705: //literal
{
feature->put(name, tr_->transcode(buf));
break;
}
case 1042: //bpchar
{
std::string str = mapnik::util::trim_copy(buf);
feature->put(name, tr_->transcode(str.c_str()));
break;
}
case 1700: //numeric
{
double val;
std::string str = numeric2string(buf);
if (mapnik::util::string2double(str, val))
{
feature->put(name, val);
}
break;
}
default:
{
MAPNIK_LOG_WARN(pgraster) << "pgraster_featureset: Unknown type_oid=" << oid;
break;
}
}
}
}
return feature;
}
return feature_ptr();
}
pgraster_featureset::~pgraster_featureset()
{
rs_->close();
}
std::string numeric2string(const char* buf)
{
boost::int16_t ndigits = int2net(buf);
boost::int16_t weight = int2net(buf+2);
boost::int16_t sign = int2net(buf+4);
boost::int16_t dscale = int2net(buf+6);
boost::scoped_array<boost::int16_t> digits(new boost::int16_t[ndigits]);
for (int n=0; n < ndigits ;++n)
{
digits[n] = int2net(buf+8+n*2);
}
std::ostringstream ss;
if (sign == 0x4000) ss << "-";
int i = std::max(weight,boost::int16_t(0));
int d = 0;
// Each numeric "digit" is actually a value between 0000 and 9999 stored in a 16 bit field.
// For example, the number 1234567809990001 is stored as four digits: [1234] [5678] [999] [1].
// Note that the last two digits show that the leading 0's are lost when the number is split.
// We must be careful to re-insert these 0's when building the string.
while ( i >= 0)
{
if (i <= weight && d < ndigits)
{
// All digits after the first must be padded to make the field 4 characters long
if (d != 0)
{
#ifdef _WINDOWS
int dig = digits[d];
if (dig < 10)
{
ss << "000"; // 0000 - 0009
}
else if (dig < 100)
{
ss << "00"; // 0010 - 0099
}
else
{
ss << "0"; // 0100 - 0999;
}
#else
switch(digits[d])
{
case 0 ... 9:
ss << "000"; // 0000 - 0009
break;
case 10 ... 99:
ss << "00"; // 0010 - 0099
break;
case 100 ... 999:
ss << "0"; // 0100 - 0999
break;
}
#endif
}
ss << digits[d++];
}
else
{
if (d == 0)
ss << "0";
else
ss << "0000";
}
i--;
}
if (dscale > 0)
{
ss << '.';
// dscale counts the number of decimal digits following the point, not the numeric digits
while (dscale > 0)
{
int value;
if (i <= weight && d < ndigits)
value = digits[d++];
else
value = 0;
// Output up to 4 decimal digits for this value
if (dscale > 0) {
ss << (value / 1000);
value %= 1000;
dscale--;
}
if (dscale > 0) {
ss << (value / 100);
value %= 100;
dscale--;
}
if (dscale > 0) {
ss << (value / 10);
value %= 10;
dscale--;
}
if (dscale > 0) {
ss << value;
dscale--;
}
i--;
}
}
return ss.str();
}

View file

@ -0,0 +1,71 @@
/*****************************************************************************
*
* This file is part of Mapnik (c++ mapping toolkit)
*
* Copyright (C) 2014 Artem Pavlenko
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
*****************************************************************************
*
* Initially developed by Sandro Santilli <strk@keybit.net>
*
*****************************************************************************/
#ifndef PGRASTER_FEATURESET_HPP
#define PGRASTER_FEATURESET_HPP
// mapnik
#include <mapnik/box2d.hpp>
#include <mapnik/datasource.hpp>
#include <mapnik/feature.hpp>
#include <mapnik/unicode.hpp>
// boost
#include <boost/scoped_ptr.hpp>
using mapnik::Featureset;
using mapnik::box2d;
using mapnik::feature_ptr;
using mapnik::raster_ptr;
using mapnik::transcoder;
using mapnik::context_ptr;
class IResultSet;
class pgraster_featureset : public mapnik::Featureset
{
public:
/// @param bandno band number (1-based). 0 (default) reads all bands.
/// Anything else forces interpretation of colors off
/// (values copied verbatim)
pgraster_featureset(boost::shared_ptr<IResultSet> const& rs,
context_ptr const& ctx,
std::string const& encoding,
bool key_field = false,
int bandno = 0);
feature_ptr next();
~pgraster_featureset();
private:
boost::shared_ptr<IResultSet> rs_;
context_ptr ctx_;
boost::scoped_ptr<mapnik::transcoder> tr_;
mapnik::value_integer feature_id_;
bool key_field_;
int band_;
};
#endif // PGRASTER_FEATURESET_HPP

View file

@ -0,0 +1,497 @@
/*****************************************************************************
*
* This file is part of Mapnik (c++ mapping toolkit)
*
* Copyright (C) 2014 Artem Pavlenko
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
*****************************************************************************
*
* Initially developed by Sandro Santilli <strk@keybit.net>
*
*****************************************************************************/
#include "pgraster_wkb_reader.hpp"
// mapnik
#include <mapnik/datasource.hpp> // for datasource_exception
#include <mapnik/global.hpp>
#include <mapnik/debug.hpp>
#include <mapnik/ctrans.hpp>
#include <mapnik/raster.hpp>
#include <mapnik/image_data.hpp>
#include <mapnik/util/conversions.hpp>
#include <mapnik/util/trim.hpp>
#include <mapnik/box2d.hpp> // for box2d
// boost
#include <boost/cstdint.hpp> // for boost::int8_t
#include <boost/make_shared.hpp>
#include <boost/bind.hpp>
namespace {
uint8_t
read_uint8(const uint8_t** from) {
return *(*from)++;
}
uint16_t
read_uint16(const boost::uint8_t** from, boost::uint8_t littleEndian) {
uint16_t ret = 0;
if (littleEndian) {
ret = (*from)[0] |
(*from)[1] << 8;
} else {
/* big endian */
ret = (*from)[0] << 8 |
(*from)[1];
}
*from += 2;
return ret;
}
int16_t
read_int16(const uint8_t** from, uint8_t littleEndian) {
assert(NULL != from);
return read_uint16(from, littleEndian);
}
double
read_float64(const boost::uint8_t** from, boost::uint8_t littleEndian) {
union {
double d;
uint64_t i;
} ret;
if (littleEndian) {
ret.i = (uint64_t) ((*from)[0] & 0xff) |
(uint64_t) ((*from)[1] & 0xff) << 8 |
(uint64_t) ((*from)[2] & 0xff) << 16 |
(uint64_t) ((*from)[3] & 0xff) << 24 |
(uint64_t) ((*from)[4] & 0xff) << 32 |
(uint64_t) ((*from)[5] & 0xff) << 40 |
(uint64_t) ((*from)[6] & 0xff) << 48 |
(uint64_t) ((*from)[7] & 0xff) << 56;
} else {
/* big endian */
ret.i = (uint64_t) ((*from)[7] & 0xff) |
(uint64_t) ((*from)[6] & 0xff) << 8 |
(uint64_t) ((*from)[5] & 0xff) << 16 |
(uint64_t) ((*from)[4] & 0xff) << 24 |
(uint64_t) ((*from)[3] & 0xff) << 32 |
(uint64_t) ((*from)[2] & 0xff) << 40 |
(uint64_t) ((*from)[1] & 0xff) << 48 |
(uint64_t) ((*from)[0] & 0xff) << 56;
}
*from += 8;
return ret.d;
}
uint32_t
read_uint32(const boost::uint8_t** from, boost::uint8_t littleEndian) {
uint32_t ret = 0;
if (littleEndian) {
ret = (uint32_t) ((*from)[0] & 0xff) |
(uint32_t) ((*from)[1] & 0xff) << 8 |
(uint32_t) ((*from)[2] & 0xff) << 16 |
(uint32_t) ((*from)[3] & 0xff) << 24;
} else {
/* big endian */
ret = (uint32_t) ((*from)[3] & 0xff) |
(uint32_t) ((*from)[2] & 0xff) << 8 |
(uint32_t) ((*from)[1] & 0xff) << 16 |
(uint32_t) ((*from)[0] & 0xff) << 24;
}
*from += 4;
return ret;
}
int32_t
read_int32(const boost::uint8_t** from, boost::uint8_t littleEndian) {
return read_uint32(from, littleEndian);
}
float
read_float32(const uint8_t** from, uint8_t littleEndian) {
union {
float f;
uint32_t i;
} ret;
ret.i = read_uint32(from, littleEndian);
return ret.f;
}
typedef enum {
PT_1BB=0, /* 1-bit boolean */
PT_2BUI=1, /* 2-bit unsigned integer */
PT_4BUI=2, /* 4-bit unsigned integer */
PT_8BSI=3, /* 8-bit signed integer */
PT_8BUI=4, /* 8-bit unsigned integer */
PT_16BSI=5, /* 16-bit signed integer */
PT_16BUI=6, /* 16-bit unsigned integer */
PT_32BSI=7, /* 32-bit signed integer */
PT_32BUI=8, /* 32-bit unsigned integer */
PT_32BF=10, /* 32-bit float */
PT_64BF=11, /* 64-bit float */
PT_END=13
} rt_pixtype;
#define BANDTYPE_FLAGS_MASK 0xF0
#define BANDTYPE_PIXTYPE_MASK 0x0F
#define BANDTYPE_FLAG_OFFDB (1<<7)
#define BANDTYPE_FLAG_HASNODATA (1<<6)
#define BANDTYPE_FLAG_ISNODATA (1<<5)
#define BANDTYPE_FLAG_RESERVED3 (1<<4)
#define BANDTYPE_PIXTYPE_MASK 0x0F
#define BANDTYPE_PIXTYPE(x) ((x)&BANDTYPE_PIXTYPE_MASK)
#define BANDTYPE_IS_OFFDB(x) ((x)&BANDTYPE_FLAG_OFFDB)
#define BANDTYPE_HAS_NODATA(x) ((x)&BANDTYPE_FLAG_HASNODATA)
#define BANDTYPE_IS_NODATA(x) ((x)&BANDTYPE_FLAG_ISNODATA)
}
using mapnik::box2d;
template<typename T>
void read_data_band(mapnik::raster_ptr raster,
uint16_t width, uint16_t height,
bool hasnodata, T reader)
{
mapnik::image_data_32 & image = raster->data_;
// Start with plain white (ABGR or RGBA depending on endiannes)
// TODO: set to transparent instead?
image.set(0xffffffff);
raster->premultiplied_alpha_ = true;
float* data = (float*)image.getBytes();
double val;
val = reader(); // nodata value, need to read anyway
if ( hasnodata ) raster->set_nodata(val);
for (int y=0; y<height; ++y) {
for (int x=0; x<width; ++x) {
val = reader();
int off = y * width + x;
data[off] = val;
}
}
}
void
pgraster_wkb_reader::read_indexed(mapnik::raster_ptr raster)
{
mapnik::image_data_32 & image = raster->data_;
// Start with all zeroes
image.set(0);
uint8_t type = read_uint8(&ptr_);
int pixtype = BANDTYPE_PIXTYPE(type);
int offline = BANDTYPE_IS_OFFDB(type) ? 1 : 0;
int hasnodata = BANDTYPE_HAS_NODATA(type) ? 1 : 0;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: band type:"
<< pixtype << " offline:" << offline
<< " hasnodata:" << hasnodata;
if ( offline ) {
MAPNIK_LOG_WARN(pgraster) << "pgraster_wkb_reader: offline band "
" unsupported";
return;
}
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: reading " << height_ << "x" << width_ << " pixels";
switch (pixtype) {
case PT_1BB:
case PT_2BUI:
case PT_4BUI:
// all <8BPP values are wrote in full bytes anyway
case PT_8BSI:
// mapnik does not support signed anyway
case PT_8BUI:
read_data_band(raster, width_, height_, hasnodata,
boost::bind(read_uint8, &ptr_));
break;
case PT_16BSI:
// mapnik does not support signed anyway
case PT_16BUI:
read_data_band(raster, width_, height_, hasnodata,
boost::bind(read_uint16, &ptr_, endian_));
break;
case PT_32BSI:
// mapnik does not support signed anyway
case PT_32BUI:
read_data_band(raster, width_, height_, hasnodata,
boost::bind(read_uint32, &ptr_, endian_));
break;
case PT_32BF:
read_data_band(raster, width_, height_, hasnodata,
boost::bind(read_float32, &ptr_, endian_));
break;
case PT_64BF:
read_data_band(raster, width_, height_, hasnodata,
boost::bind(read_float64, &ptr_, endian_));
break;
default:
std::ostringstream err;
err << "pgraster_wkb_reader: data band type " << pixtype << " unsupported";
// TODO: accept policy to decide on throw-or-skip ?
//MAPNIK_LOG_WARN(pgraster) << err.str();
throw mapnik::datasource_exception(err.str());
}
}
template<typename T>
void read_grayscale_band(mapnik::raster_ptr raster,
uint16_t width, uint16_t height,
bool hasnodata, T reader)
{
mapnik::image_data_32 & image = raster->data_;
// Start with plain white (ABGR or RGBA depending on endiannes)
// TODO: set to transparent instead?
image.set(0xffffffff);
raster->premultiplied_alpha_ = true;
int val;
uint8_t * data = image.getBytes();
int ps = 4; // sizeof(image_data::pixel_type)
int off;
val = reader(); // nodata value, need to read anyway
if ( hasnodata ) raster->set_nodata(val);
for (int y=0; y<height; ++y) {
for (int x=0; x<width; ++x) {
val = reader();
off = y * width * ps + x * ps;
// Pixel space is RGBA
data[off+0] = val;
data[off+1] = val;
data[off+2] = val;
}
}
}
void
pgraster_wkb_reader::read_grayscale(mapnik::raster_ptr raster)
{
uint8_t type = read_uint8(&ptr_);
int pixtype = BANDTYPE_PIXTYPE(type);
int offline = BANDTYPE_IS_OFFDB(type) ? 1 : 0;
int hasnodata = BANDTYPE_HAS_NODATA(type) ? 1 : 0;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: band type:"
<< pixtype << " offline:" << offline
<< " hasnodata:" << hasnodata;
if ( offline ) {
MAPNIK_LOG_WARN(pgraster) << "pgraster_wkb_reader: offline band "
" unsupported";
return;
}
switch (pixtype) {
case PT_1BB:
case PT_2BUI:
case PT_4BUI:
// all <8BPP values are wrote in full bytes anyway
case PT_8BSI:
// mapnik does not support signed anyway
case PT_8BUI:
read_grayscale_band(raster, width_, height_, hasnodata,
boost::bind(read_uint8, &ptr_));
break;
case PT_16BSI:
// mapnik does not support signed anyway
case PT_16BUI:
read_grayscale_band(raster, width_, height_, hasnodata,
boost::bind(read_uint16, &ptr_, endian_));
break;
case PT_32BSI:
// mapnik does not support signed anyway
case PT_32BUI:
read_grayscale_band(raster, width_, height_, hasnodata,
boost::bind(read_uint32, &ptr_, endian_));
break;
default:
std::ostringstream err;
err << "pgraster_wkb_reader: grayscale band type "
<< pixtype << " unsupported";
//MAPNIK_LOG_WARN(pgraster) << err.str();
throw mapnik::datasource_exception(err.str());
}
}
void
pgraster_wkb_reader::read_rgba(mapnik::raster_ptr raster)
{
mapnik::image_data_32 & image = raster->data_;
// Start with plain white (ABGR or RGBA depending on endiannes)
image.set(0xffffffff);
//raster->set_nodata(0xffffffff);
raster->premultiplied_alpha_ = true;
uint8_t nodataval;
for (int bn=0; bn<numBands_; ++bn) {
uint8_t type = read_uint8(&ptr_);
int pixtype = BANDTYPE_PIXTYPE(type);
int offline = BANDTYPE_IS_OFFDB(type) ? 1 : 0;
int hasnodata = BANDTYPE_HAS_NODATA(type) ? 1 : 0;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: band " << bn
<< " type:" << pixtype << " offline:" << offline
<< " hasnodata:" << hasnodata;
if ( offline ) {
MAPNIK_LOG_WARN(pgraster) << "pgraster_wkb_reader: offline band "
<< bn << " unsupported";
continue;
}
if ( pixtype > PT_8BUI || pixtype < PT_8BSI ) {
MAPNIK_LOG_WARN(pgraster) << "pgraster_wkb_reader: band " << bn
<< " type " << type << " unsupported";
continue;
}
uint8_t tmp = read_uint8(&ptr_);
if ( ! bn ) nodataval = tmp;
else if ( tmp != nodataval ) {
MAPNIK_LOG_WARN(pgraster) << "pgraster_wkb_reader: band " << bn
<< " nodataval " << tmp << " != band 0 nodataval " << nodataval;
}
int ps = 4; // sizeof(image_data::pixel_type)
uint8_t * image_data = image.getBytes();
for (int y=0; y<height_; ++y) {
for (int x=0; x<width_; ++x) {
uint8_t val = read_uint8(&ptr_);
// y * width_ * ps is the row (ps is pixel size)
// x * ps is the column
int off = y * width_ * ps + x * ps;
// Pixel space is RGBA
image_data[off+bn] = val;
}
}
}
}
mapnik::raster_ptr
pgraster_wkb_reader::get_raster() {
/* Read endianness */
endian_ = *ptr_;
ptr_ += 1;
/* Read version of protocol */
uint16_t version = read_uint16(&ptr_, endian_);
if (version != 0) {
MAPNIK_LOG_WARN(pgraster) << "pgraster_wkb_reader: WKB version "
<< version << " unsupported";
return mapnik::raster_ptr();
}
numBands_ = read_uint16(&ptr_, endian_);
double scaleX = read_float64(&ptr_, endian_);
double scaleY = read_float64(&ptr_, endian_);
double ipX = read_float64(&ptr_, endian_);
double ipY = read_float64(&ptr_, endian_);
double skewX = read_float64(&ptr_, endian_);
double skewY = read_float64(&ptr_, endian_);
int32_t srid = read_int32(&ptr_, endian_);
width_ = read_uint16(&ptr_, endian_);
height_ = read_uint16(&ptr_, endian_);
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: numBands=" << numBands_;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: scaleX=" << scaleX;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: scaleY=" << scaleY;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: ipX=" << ipX;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: ipY=" << ipY;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: skewX=" << skewX;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: skewY=" << skewY;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: srid=" << srid;
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: size="
<< width_ << "x" << height_;
// this is for color interpretation
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: bandno=" << bandno_;
if ( skewX || skewY ) {
MAPNIK_LOG_WARN(pgraster) << "pgraster_wkb_reader: raster rotation is not supported";
return mapnik::raster_ptr();
}
box2d<double> ext(ipX,ipY,ipX+(width_*scaleX),ipY+(height_*scaleY));
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: Raster extent=" << ext;
mapnik::raster_ptr raster = boost::make_shared<mapnik::raster>(ext, width_, height_);
if ( bandno_ ) {
if ( bandno_ != 1 ) {
MAPNIK_LOG_WARN(pgraster) << "pgraster_wkb_reader: "
"reading bands other than 1st as indexed is unsupported";
return mapnik::raster_ptr();
}
MAPNIK_LOG_DEBUG(pgraster) << "pgraster_wkb_reader: requested band " << bandno_;
read_indexed(raster);
}
else {
switch (numBands_) {
case 1:
read_grayscale(raster);
break;
case 3:
case 4:
read_rgba(raster);
break;
default:
std::ostringstream err;
err << "pgraster_wkb_reader: raster with "
<< numBands_
<< " bands is not supported, specify a band number";
//MAPNIK_LOG_WARN(pgraster) << err.str();
throw mapnik::datasource_exception(err.str());
return mapnik::raster_ptr();
}
}
return raster;
}

View file

@ -0,0 +1,87 @@
/*****************************************************************************
*
* This file is part of Mapnik (c++ mapping toolkit)
*
* Copyright (C) 2014 Artem Pavlenko
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
*****************************************************************************
*
* Initially developed by Sandro Santilli <strk@keybit.net>
*
*****************************************************************************/
#ifndef PGRASTER_WKB_READER_HPP
#define PGRASTER_WKB_READER_HPP
// mapnik
#include <mapnik/feature.hpp> // for raster_ptr
// boost
#include <boost/cstdint.hpp> // for boost::uint8_t
enum pgraster_color_interp {
// Automatic color interpretation:
// uses grayscale for single band, rgb for 3 bands
// rgba for 4 bands
pgr_auto,
// Grayscale interpretation
pgr_grayscale,
pgr_indexed,
pgr_rgb,
pgr_rgba
};
class pgraster_wkb_reader
{
public:
pgraster_wkb_reader(const uint8_t* wkb, int size, int bnd=0)
: wkbsize_(size), wkb_(wkb), wkbend_(wkb+size), ptr_(wkb), bandno_(bnd)
{}
mapnik::raster_ptr get_raster();
/// @param bnd band number. If 0 (default) it'll try to read all bands
/// with automatic color interpretation (rgb for 3 bands,
/// rgba for 4 bands, grayscale for 1 band).
/// Any other value results in pixel
/// values being copied verbatim into the returned raster
/// for interpretation by the caller.
static mapnik::raster_ptr read(const uint8_t* wkb, int size, int bnd=0)
{
pgraster_wkb_reader reader(wkb,size,bnd);
return reader.get_raster();
}
private:
void read_indexed(mapnik::raster_ptr raster);
void read_grayscale(mapnik::raster_ptr raster);
void read_rgba(mapnik::raster_ptr raster);
int wkbsize_;
const uint8_t* wkb_;
const uint8_t* wkbend_;
const uint8_t* ptr_;
uint8_t endian_;
int bandno_;
uint16_t numBands_;
uint16_t width_;
uint16_t height_;
};
#endif // PGRASTER_WKB_READER_HPP

View file

@ -0,0 +1,747 @@
#!/usr/bin/env python
from nose.tools import *
import atexit
import cProfile, pstats, io
import time
from utilities import execution_path, run_all
from subprocess import Popen, PIPE
import os, mapnik
from Queue import Queue
import threading
import sys
import re
from binascii import hexlify, unhexlify
MAPNIK_TEST_DBNAME = 'mapnik-tmp-pgraster-test-db'
POSTGIS_TEMPLATE_DBNAME = 'template_postgis'
def setup():
# All of the paths used are relative, if we run the tests
# from another directory we need to chdir()
os.chdir(execution_path('.'))
def call(cmd,silent=False):
proc = Popen(cmd, shell=True, stdout=PIPE, stderr=PIPE)
stdout, stderr = proc.communicate()
if stderr and not silent:
print stderr.strip()
if proc.returncode:
raise RuntimeError(stderr.strip())
def psql_can_connect():
"""Test ability to connect to a postgis template db with no options.
Basically, to run these tests your user must have full read
access over unix sockets without supplying a password. This
keeps these tests simple and focused on postgis not on postgres
auth issues.
"""
try:
call('psql %s -c "select postgis_version()"' % POSTGIS_TEMPLATE_DBNAME)
return True
except RuntimeError, e:
print 'Notice: skipping postgis tests (connection)'
return False
def psql_run(cmd):
cmd = 'psql --set ON_ERROR_STOP=1 %s -c "%s"' % \
(MAPNIK_TEST_DBNAME, cmd.replace('"', '\\"'))
print 'DEBUG: running ' + cmd
call(cmd)
def raster2pgsql_on_path():
"""Test for presence of raster2pgsql on the user path.
We require this program to load test data into a temporarily database.
"""
try:
call('raster2pgsql')
return True
except RuntimeError, e:
print 'Notice: skipping postgis tests (raster2pgsql)'
return False
def createdb_and_dropdb_on_path():
"""Test for presence of dropdb/createdb on user path.
We require these programs to setup and teardown the testing db.
"""
try:
call('createdb --help')
call('dropdb --help')
return True
except RuntimeError, e:
print 'Notice: skipping postgis tests (createdb/dropdb)'
return False
def postgis_setup():
call('dropdb %s' % MAPNIK_TEST_DBNAME,silent=True)
call('createdb -T %s %s' % (POSTGIS_TEMPLATE_DBNAME,MAPNIK_TEST_DBNAME),silent=False)
def postgis_takedown():
pass
# fails as the db is in use: https://github.com/mapnik/mapnik/issues/960
#call('dropdb %s' % MAPNIK_TEST_DBNAME)
def import_raster(filename, tabname, tilesize, constraint, overview):
print 'tile: ' + tilesize + ' constraints: ' + str(constraint) \
+ ' overviews: ' + overview
cmd = 'raster2pgsql -Y -I -q'
if constraint:
cmd += ' -C'
if tilesize:
cmd += ' -t ' + tilesize
if overview:
cmd += ' -l ' + overview
cmd += ' %s %s | psql --set ON_ERROR_STOP=1 -q %s' % (filename,tabname,MAPNIK_TEST_DBNAME)
print 'Import call: ' + cmd
call(cmd)
def drop_imported(tabname, overview):
psql_run('DROP TABLE IF EXISTS "' + tabname + '";')
if overview:
for of in overview.split(','):
psql_run('DROP TABLE IF EXISTS "o_' + of + '_' + tabname + '";')
if 'pgraster' in mapnik.DatasourceCache.plugin_names() \
and createdb_and_dropdb_on_path() \
and psql_can_connect() \
and raster2pgsql_on_path():
# initialize test database
postgis_setup()
# dataraster.tif, 2283x1913 int16 single-band
def _test_dataraster_16bsi_rendering(lbl, overview, rescale, clip):
if rescale:
lbl += ' Sc'
if clip:
lbl += ' Cl'
ds = mapnik.PgRaster(dbname=MAPNIK_TEST_DBNAME,table='"dataRaster"',
band=1,use_overviews=1 if overview else 0,
prescale_rasters=rescale,clip_rasters=clip)
fs = ds.featureset()
feature = fs.next()
eq_(feature['rid'],1)
lyr = mapnik.Layer('dataraster_16bsi')
lyr.datasource = ds
expenv = mapnik.Box2d(-14637, 3903178, 1126863, 4859678)
env = lyr.envelope()
# As the input size is a prime number both horizontally
# and vertically, we expect the extent of the overview
# tables to be a pixel wider than the original, whereas
# the pixel size in geographical units depends on the
# overview factor. So we start with the original pixel size
# as base scale and multiply by the overview factor.
# NOTE: the overview table extent only grows north and east
pixsize = 500 # see gdalinfo dataraster.tif
tol = pixsize * max(overview.split(',')) if overview else 0
assert_almost_equal(env.minx, expenv.minx)
assert_almost_equal(env.miny, expenv.miny, delta=tol)
assert_almost_equal(env.maxx, expenv.maxx, delta=tol)
assert_almost_equal(env.maxy, expenv.maxy)
mm = mapnik.Map(256, 256)
style = mapnik.Style()
col = mapnik.RasterColorizer();
col.default_mode = mapnik.COLORIZER_DISCRETE;
col.add_stop(0, mapnik.Color(0x40,0x40,0x40,255));
col.add_stop(10, mapnik.Color(0x80,0x80,0x80,255));
col.add_stop(20, mapnik.Color(0xa0,0xa0,0xa0,255));
sym = mapnik.RasterSymbolizer()
sym.colorizer = col
rule = mapnik.Rule()
rule.symbols.append(sym)
style.rules.append(rule)
mm.append_style('foo', style)
lyr.styles.append('foo')
mm.layers.append(lyr)
mm.zoom_to_box(expenv)
im = mapnik.Image(mm.width, mm.height)
t0 = time.time() # we want wall time to include IO waits
mapnik.render(mm, im)
lap = time.time() - t0
print 'T ' + str(lap) + ' -- ' + lbl + ' E:full'
# no data
eq_(im.view(1,1,1,1).tostring(), '\x00\x00\x00\x00')
eq_(im.view(255,255,1,1).tostring(), '\x00\x00\x00\x00')
eq_(im.view(195,116,1,1).tostring(), '\x00\x00\x00\x00')
# A0A0A0
eq_(im.view(100,120,1,1).tostring(), '\xa0\xa0\xa0\xff')
eq_(im.view( 75, 80,1,1).tostring(), '\xa0\xa0\xa0\xff')
# 808080
eq_(im.view( 74,170,1,1).tostring(), '\x80\x80\x80\xff')
eq_(im.view( 30, 50,1,1).tostring(), '\x80\x80\x80\xff')
# 404040
eq_(im.view(190, 70,1,1).tostring(), '\x40\x40\x40\xff')
eq_(im.view(140,170,1,1).tostring(), '\x40\x40\x40\xff')
# Now zoom over a portion of the env (1/10)
newenv = mapnik.Box2d(273663,4024478,330738,4072303)
mm.zoom_to_box(newenv)
t0 = time.time() # we want wall time to include IO waits
mapnik.render(mm, im)
lap = time.time() - t0
print 'T ' + str(lap) + ' -- ' + lbl + ' E:1/10'
# nodata
eq_(hexlify(im.view(255,255,1,1).tostring()), '00000000')
eq_(hexlify(im.view(200,254,1,1).tostring()), '00000000')
# A0A0A0
eq_(hexlify(im.view(90,232,1,1).tostring()), 'a0a0a0ff')
eq_(hexlify(im.view(96,245,1,1).tostring()), 'a0a0a0ff')
# 808080
eq_(hexlify(im.view(1,1,1,1).tostring()), '808080ff')
eq_(hexlify(im.view(128,128,1,1).tostring()), '808080ff')
# 404040
eq_(hexlify(im.view(255, 0,1,1).tostring()), '404040ff')
def _test_dataraster_16bsi(lbl, tilesize, constraint, overview):
rf = os.path.join(execution_path('.'),'../data/raster/dataraster.tif')
import_raster(rf, 'dataRaster', tilesize, constraint, overview)
if constraint:
lbl += ' C'
if tilesize:
lbl += ' T:' + tilesize
if overview:
lbl += ' O:' + overview
for prescale in [0,1]:
for clip in [0,1]:
_test_dataraster_16bsi_rendering(lbl, overview, prescale, clip)
drop_imported('dataRaster', overview)
def test_dataraster_16bsi():
for tilesize in ['','256x256']:
for constraint in [0,1]:
for overview in ['','4','2,16']:
_test_dataraster_16bsi('data_16bsi', tilesize, constraint, overview)
# river.tiff, RGBA 8BUI
def _test_rgba_8bui_rendering(lbl, overview, rescale, clip):
if rescale:
lbl += ' Sc'
if clip:
lbl += ' Cl'
ds = mapnik.PgRaster(dbname=MAPNIK_TEST_DBNAME,table='(select * from "River") foo',
use_overviews=1 if overview else 0,
prescale_rasters=rescale,clip_rasters=clip)
fs = ds.featureset()
feature = fs.next()
eq_(feature['rid'],1)
lyr = mapnik.Layer('rgba_8bui')
lyr.datasource = ds
expenv = mapnik.Box2d(0, -210, 256, 0)
env = lyr.envelope()
# As the input size is a prime number both horizontally
# and vertically, we expect the extent of the overview
# tables to be a pixel wider than the original, whereas
# the pixel size in geographical units depends on the
# overview factor. So we start with the original pixel size
# as base scale and multiply by the overview factor.
# NOTE: the overview table extent only grows north and east
pixsize = 1 # see gdalinfo river.tif
tol = pixsize * max(overview.split(',')) if overview else 0
assert_almost_equal(env.minx, expenv.minx)
assert_almost_equal(env.miny, expenv.miny, delta=tol)
assert_almost_equal(env.maxx, expenv.maxx, delta=tol)
assert_almost_equal(env.maxy, expenv.maxy)
mm = mapnik.Map(256, 256)
style = mapnik.Style()
sym = mapnik.RasterSymbolizer()
rule = mapnik.Rule()
rule.symbols.append(sym)
style.rules.append(rule)
mm.append_style('foo', style)
lyr.styles.append('foo')
mm.layers.append(lyr)
mm.zoom_to_box(expenv)
im = mapnik.Image(mm.width, mm.height)
t0 = time.time() # we want wall time to include IO waits
mapnik.render(mm, im)
lap = time.time() - t0
print 'T ' + str(lap) + ' -- ' + lbl + ' E:full'
#im.save('/tmp/xfull.png') # for debugging
# no data
eq_(hexlify(im.view(3,3,1,1).tostring()), '00000000')
eq_(hexlify(im.view(250,250,1,1).tostring()), '00000000')
# full opaque river color
eq_(hexlify(im.view(175,118,1,1).tostring()), 'b9d8f8ff')
# half-transparent pixel
pxstr = hexlify(im.view(122,138,1,1).tostring())
apat = ".*(..)$"
match = re.match(apat, pxstr)
assert match, 'pixel ' + pxstr + ' does not match pattern ' + apat
alpha = match.group(1)
assert alpha != 'ff' and alpha != '00', \
'unexpected full transparent/opaque pixel: ' + alpha
# Now zoom over a portion of the env (1/10)
newenv = mapnik.Box2d(166,-105,191,-77)
mm.zoom_to_box(newenv)
t0 = time.time() # we want wall time to include IO waits
mapnik.render(mm, im)
lap = time.time() - t0
print 'T ' + str(lap) + ' -- ' + lbl + ' E:1/10'
#im.save('/tmp/xtenth.png') # for debugging
# no data
eq_(hexlify(im.view(255,255,1,1).tostring()), '00000000')
eq_(hexlify(im.view(200,40,1,1).tostring()), '00000000')
# full opaque river color
eq_(hexlify(im.view(100,168,1,1).tostring()), 'b9d8f8ff')
# half-transparent pixel
pxstr = hexlify(im.view(122,138,1,1).tostring())
apat = ".*(..)$"
match = re.match(apat, pxstr)
assert match, 'pixel ' + pxstr + ' does not match pattern ' + apat
alpha = match.group(1)
assert alpha != 'ff' and alpha != '00', \
'unexpected full transparent/opaque pixel: ' + alpha
def _test_rgba_8bui(lbl, tilesize, constraint, overview):
rf = os.path.join(execution_path('.'),'../data/raster/river.tiff')
import_raster(rf, 'River', tilesize, constraint, overview)
if constraint:
lbl += ' C'
if tilesize:
lbl += ' T:' + tilesize
if overview:
lbl += ' O:' + overview
for prescale in [0,1]:
for clip in [0,1]:
_test_rgba_8bui_rendering(lbl, overview, prescale, clip)
drop_imported('River', overview)
def test_rgba_8bui():
for tilesize in ['','16x16']:
for constraint in [0,1]:
for overview in ['2']:
_test_rgba_8bui('rgba_8bui', tilesize, constraint, overview)
# nodata-edge.tif, RGB 8BUI
def _test_rgb_8bui_rendering(lbl, tnam, overview, rescale, clip):
if rescale:
lbl += ' Sc'
if clip:
lbl += ' Cl'
ds = mapnik.PgRaster(dbname=MAPNIK_TEST_DBNAME,table=tnam,
use_overviews=1 if overview else 0,
prescale_rasters=rescale,clip_rasters=clip)
fs = ds.featureset()
feature = fs.next()
eq_(feature['rid'],1)
lyr = mapnik.Layer('rgba_8bui')
lyr.datasource = ds
expenv = mapnik.Box2d(-12329035.7652168,4508650.39854396, \
-12328653.0279471,4508957.34625536)
env = lyr.envelope()
# As the input size is a prime number both horizontally
# and vertically, we expect the extent of the overview
# tables to be a pixel wider than the original, whereas
# the pixel size in geographical units depends on the
# overview factor. So we start with the original pixel size
# as base scale and multiply by the overview factor.
# NOTE: the overview table extent only grows north and east
pixsize = 2 # see gdalinfo nodata-edge.tif
tol = pixsize * max(overview.split(',')) if overview else 0
assert_almost_equal(env.minx, expenv.minx, places=0)
assert_almost_equal(env.miny, expenv.miny, delta=tol)
assert_almost_equal(env.maxx, expenv.maxx, delta=tol)
assert_almost_equal(env.maxy, expenv.maxy, places=0)
mm = mapnik.Map(256, 256)
style = mapnik.Style()
sym = mapnik.RasterSymbolizer()
rule = mapnik.Rule()
rule.symbols.append(sym)
style.rules.append(rule)
mm.append_style('foo', style)
lyr.styles.append('foo')
mm.layers.append(lyr)
mm.zoom_to_box(expenv)
im = mapnik.Image(mm.width, mm.height)
t0 = time.time() # we want wall time to include IO waits
mapnik.render(mm, im)
lap = time.time() - t0
print 'T ' + str(lap) + ' -- ' + lbl + ' E:full'
#im.save('/tmp/xfull.png') # for debugging
# no data
eq_(hexlify(im.view(3,16,1,1).tostring()), '00000000')
eq_(hexlify(im.view(128,16,1,1).tostring()), '00000000')
eq_(hexlify(im.view(250,16,1,1).tostring()), '00000000')
eq_(hexlify(im.view(3,240,1,1).tostring()), '00000000')
eq_(hexlify(im.view(128,240,1,1).tostring()), '00000000')
eq_(hexlify(im.view(250,240,1,1).tostring()), '00000000')
# dark brown
eq_(hexlify(im.view(174,39,1,1).tostring()), 'c3a698ff')
# dark gray
eq_(hexlify(im.view(195,132,1,1).tostring()), '575f62ff')
# Now zoom over a portion of the env (1/10)
newenv = mapnik.Box2d(-12329035.7652168, 4508926.651484220, \
-12328997.49148983,4508957.34625536)
mm.zoom_to_box(newenv)
t0 = time.time() # we want wall time to include IO waits
mapnik.render(mm, im)
lap = time.time() - t0
print 'T ' + str(lap) + ' -- ' + lbl + ' E:1/10'
#im.save('/tmp/xtenth.png') # for debugging
# no data
eq_(hexlify(im.view(3,16,1,1).tostring()), '00000000')
eq_(hexlify(im.view(128,16,1,1).tostring()), '00000000')
eq_(hexlify(im.view(250,16,1,1).tostring()), '00000000')
# black
eq_(hexlify(im.view(3,42,1,1).tostring()), '000000ff')
eq_(hexlify(im.view(3,134,1,1).tostring()), '000000ff')
eq_(hexlify(im.view(3,244,1,1).tostring()), '000000ff')
# gray
eq_(hexlify(im.view(135,157,1,1).tostring()), '4e555bff')
# brown
eq_(hexlify(im.view(195,223,1,1).tostring()), 'f2cdbaff')
def _test_rgb_8bui(lbl, tilesize, constraint, overview):
rf = os.path.join(execution_path('.'),'../data/raster/nodata-edge.tif')
tnam = 'nodataedge'
import_raster(rf, tnam, tilesize, constraint, overview)
if constraint:
lbl += ' C'
if tilesize:
lbl += ' T:' + tilesize
if overview:
lbl += ' O:' + overview
for prescale in [0,1]:
for clip in [0,1]:
_test_rgb_8bui_rendering(lbl, tnam, overview, prescale, clip)
#drop_imported(tnam, overview)
def test_rgb_8bui():
for tilesize in ['64x64']:
for constraint in [1]:
for overview in ['']:
_test_rgb_8bui('rgb_8bui', tilesize, constraint, overview)
def _test_grayscale_subquery(lbl,pixtype,value):
#
# 3 8 13
# +---+---+---+
# 3 | v | v | v | NOTE: writes different color
# +---+---+---+ in 13,8 and 8,13
# 8 | v | v | a |
# +---+---+---+
# 13 | v | b | v |
# +---+---+---+
#
val_a = value/3;
val_b = val_a*2;
sql = "(select 3 as i, " \
" ST_SetValues(" \
" ST_SetValues(" \
" ST_AsRaster(" \
" ST_MakeEnvelope(0,0,14,14), " \
" 1.0, -1.0, '%s', %s" \
" ), " \
" 11, 6, 4, 5, %s::float8" \
" )," \
" 6, 11, 5, 4, %s::float8" \
" ) as \"R\"" \
") as foo" % (pixtype,value, val_a, val_b)
rescale = 0
clip = 0
if rescale:
lbl += ' Sc'
if clip:
lbl += ' Cl'
ds = mapnik.PgRaster(dbname=MAPNIK_TEST_DBNAME, table=sql,
raster_field='"R"', use_overviews=1,
prescale_rasters=rescale,clip_rasters=clip)
fs = ds.featureset()
feature = fs.next()
eq_(feature['i'],3)
lyr = mapnik.Layer('grayscale_subquery')
lyr.datasource = ds
expenv = mapnik.Box2d(0,0,14,14)
env = lyr.envelope()
assert_almost_equal(env.minx, expenv.minx, places=0)
assert_almost_equal(env.miny, expenv.miny, places=0)
assert_almost_equal(env.maxx, expenv.maxx, places=0)
assert_almost_equal(env.maxy, expenv.maxy, places=0)
mm = mapnik.Map(15, 15)
style = mapnik.Style()
sym = mapnik.RasterSymbolizer()
rule = mapnik.Rule()
rule.symbols.append(sym)
style.rules.append(rule)
mm.append_style('foo', style)
lyr.styles.append('foo')
mm.layers.append(lyr)
mm.zoom_to_box(expenv)
im = mapnik.Image(mm.width, mm.height)
t0 = time.time() # we want wall time to include IO waits
mapnik.render(mm, im)
lap = time.time() - t0
print 'T ' + str(lap) + ' -- ' + lbl + ' E:full'
#im.save('/tmp/xfull.png') # for debugging
h = format(value, '02x')
hex_v = h+h+h+'ff'
h = format(val_a, '02x')
hex_a = h+h+h+'ff'
h = format(val_b, '02x')
hex_b = h+h+h+'ff'
eq_(hexlify(im.view( 3, 3,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 8, 3,1,1).tostring()), hex_v);
eq_(hexlify(im.view(13, 3,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 3, 8,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 8, 8,1,1).tostring()), hex_v);
eq_(hexlify(im.view(13, 8,1,1).tostring()), hex_a);
eq_(hexlify(im.view( 3,13,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 8,13,1,1).tostring()), hex_b);
eq_(hexlify(im.view(13,13,1,1).tostring()), hex_v);
def test_grayscale_2bui_subquery():
_test_grayscale_subquery('grayscale_2bui_subquery', '2BUI', 3)
def test_grayscale_4bui_subquery():
_test_grayscale_subquery('grayscale_4bui_subquery', '4BUI', 15)
def test_grayscale_8bui_subquery():
_test_grayscale_subquery('grayscale_8bui_subquery', '8BUI', 63)
def test_grayscale_8bsi_subquery():
# NOTE: we're using a positive integer because Mapnik
# does not support negative data values anyway
_test_grayscale_subquery('grayscale_8bsi_subquery', '8BSI', 69)
def test_grayscale_16bui_subquery():
_test_grayscale_subquery('grayscale_16bui_subquery', '16BUI', 126)
def test_grayscale_16bsi_subquery():
# NOTE: we're using a positive integer because Mapnik
# does not support negative data values anyway
_test_grayscale_subquery('grayscale_16bsi_subquery', '16BSI', 144)
def test_grayscale_32bui_subquery():
_test_grayscale_subquery('grayscale_32bui_subquery', '32BUI', 255)
def test_grayscale_32bsi_subquery():
# NOTE: we're using a positive integer because Mapnik
# does not support negative data values anyway
_test_grayscale_subquery('grayscale_32bsi_subquery', '32BSI', 129)
def _test_data_subquery(lbl, pixtype, value):
#
# 3 8 13
# +---+---+---+
# 3 | v | v | v | NOTE: writes different values
# +---+---+---+ in 13,8 and 8,13
# 8 | v | v | a |
# +---+---+---+
# 13 | v | b | v |
# +---+---+---+
#
val_a = value/3;
val_b = val_a*2;
sql = "(select 3 as i, " \
" ST_SetValues(" \
" ST_SetValues(" \
" ST_AsRaster(" \
" ST_MakeEnvelope(0,0,14,14), " \
" 1.0, -1.0, '%s', %s" \
" ), " \
" 11, 6, 5, 5, %s::float8" \
" )," \
" 6, 11, 5, 5, %s::float8" \
" ) as \"R\"" \
") as foo" % (pixtype,value, val_a, val_b)
overview = ''
rescale = 0
clip = 0
if rescale:
lbl += ' Sc'
if clip:
lbl += ' Cl'
ds = mapnik.PgRaster(dbname=MAPNIK_TEST_DBNAME, table=sql,
raster_field='R', use_overviews=0 if overview else 0,
band=1, prescale_rasters=rescale, clip_rasters=clip)
fs = ds.featureset()
feature = fs.next()
eq_(feature['i'],3)
lyr = mapnik.Layer('data_subquery')
lyr.datasource = ds
expenv = mapnik.Box2d(0,0,14,14)
env = lyr.envelope()
assert_almost_equal(env.minx, expenv.minx, places=0)
assert_almost_equal(env.miny, expenv.miny, places=0)
assert_almost_equal(env.maxx, expenv.maxx, places=0)
assert_almost_equal(env.maxy, expenv.maxy, places=0)
mm = mapnik.Map(15, 15)
style = mapnik.Style()
col = mapnik.RasterColorizer();
col.default_mode = mapnik.COLORIZER_DISCRETE;
col.add_stop(val_a, mapnik.Color(0xff,0x00,0x00,255));
col.add_stop(val_b, mapnik.Color(0x00,0xff,0x00,255));
col.add_stop(value, mapnik.Color(0x00,0x00,0xff,255));
sym = mapnik.RasterSymbolizer()
sym.colorizer = col
rule = mapnik.Rule()
rule.symbols.append(sym)
style.rules.append(rule)
mm.append_style('foo', style)
lyr.styles.append('foo')
mm.layers.append(lyr)
mm.zoom_to_box(expenv)
im = mapnik.Image(mm.width, mm.height)
t0 = time.time() # we want wall time to include IO waits
mapnik.render(mm, im)
lap = time.time() - t0
print 'T ' + str(lap) + ' -- ' + lbl + ' E:full'
#im.save('/tmp/xfull.png') # for debugging
h = format(value, '02x')
hex_v = '0000ffff'
hex_a = 'ff0000ff'
hex_b = '00ff00ff'
eq_(hexlify(im.view( 3, 3,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 8, 3,1,1).tostring()), hex_v);
eq_(hexlify(im.view(13, 3,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 3, 8,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 8, 8,1,1).tostring()), hex_v);
eq_(hexlify(im.view(13, 8,1,1).tostring()), hex_a);
eq_(hexlify(im.view( 3,13,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 8,13,1,1).tostring()), hex_b);
eq_(hexlify(im.view(13,13,1,1).tostring()), hex_v);
def test_data_2bui_subquery():
_test_data_subquery('data_2bui_subquery', '2BUI', 3)
def test_data_4bui_subquery():
_test_data_subquery('data_4bui_subquery', '4BUI', 15)
def test_data_8bui_subquery():
_test_data_subquery('data_8bui_subquery', '8BUI', 63)
def test_data_8bsi_subquery():
# NOTE: we're using a positive integer because Mapnik
# does not support negative data values anyway
_test_data_subquery('data_8bsi_subquery', '8BSI', 69)
def test_data_16bui_subquery():
_test_data_subquery('data_16bui_subquery', '16BUI', 126)
def test_data_16bsi_subquery():
# NOTE: we're using a positive integer because Mapnik
# does not support negative data values anyway
_test_data_subquery('data_16bsi_subquery', '16BSI', 135)
def test_data_32bui_subquery():
_test_data_subquery('data_32bui_subquery', '32BUI', 255)
def test_data_32bsi_subquery():
# NOTE: we're using a positive integer because Mapnik
# does not support negative data values anyway
_test_data_subquery('data_32bsi_subquery', '32BSI', 264)
def test_data_32bf_subquery():
_test_data_subquery('data_32bf_subquery', '32BF', 450)
def test_data_64bf_subquery():
_test_data_subquery('data_64bf_subquery', '64BF', 3072)
def _test_rgba_subquery(lbl, pixtype, r, g, b, a, g1, b1):
#
# 3 8 13
# +---+---+---+
# 3 | v | v | h | NOTE: writes different alpha
# +---+---+---+ in 13,8 and 8,13
# 8 | v | v | a |
# +---+---+---+
# 13 | v | b | v |
# +---+---+---+
#
sql = "(select 3 as i, " \
" ST_SetValues(" \
" ST_SetValues(" \
" ST_AddBand(" \
" ST_AddBand(" \
" ST_AddBand(" \
" ST_AsRaster(" \
" ST_MakeEnvelope(0,0,14,14), " \
" 1.0, -1.0, '%s', %s" \
" )," \
" '%s', %d::float" \
" ), " \
" '%s', %d::float" \
" ), " \
" '%s', %d::float" \
" ), " \
" 2, 11, 6, 4, 5, %s::float8" \
" )," \
" 3, 6, 11, 5, 4, %s::float8" \
" ) as r" \
") as foo" % (pixtype, r, pixtype, g, pixtype, b, pixtype, a, g1, b1)
overview = ''
rescale = 0
clip = 0
if rescale:
lbl += ' Sc'
if clip:
lbl += ' Cl'
ds = mapnik.PgRaster(dbname=MAPNIK_TEST_DBNAME, table=sql,
raster_field='r', use_overviews=0 if overview else 0,
prescale_rasters=rescale, clip_rasters=clip)
fs = ds.featureset()
feature = fs.next()
eq_(feature['i'],3)
lyr = mapnik.Layer('rgba_subquery')
lyr.datasource = ds
expenv = mapnik.Box2d(0,0,14,14)
env = lyr.envelope()
assert_almost_equal(env.minx, expenv.minx, places=0)
assert_almost_equal(env.miny, expenv.miny, places=0)
assert_almost_equal(env.maxx, expenv.maxx, places=0)
assert_almost_equal(env.maxy, expenv.maxy, places=0)
mm = mapnik.Map(15, 15)
style = mapnik.Style()
sym = mapnik.RasterSymbolizer()
rule = mapnik.Rule()
rule.symbols.append(sym)
style.rules.append(rule)
mm.append_style('foo', style)
lyr.styles.append('foo')
mm.layers.append(lyr)
mm.zoom_to_box(expenv)
im = mapnik.Image(mm.width, mm.height)
t0 = time.time() # we want wall time to include IO waits
mapnik.render(mm, im)
lap = time.time() - t0
print 'T ' + str(lap) + ' -- ' + lbl + ' E:full'
im.save('/tmp/xfull.png') # for debugging
hex_v = format(r << 24 | g << 16 | b << 8 | a, '08x')
hex_a = format(r << 24 | g1 << 16 | b << 8 | a, '08x')
hex_b = format(r << 24 | g << 16 | b1 << 8 | a, '08x')
eq_(hexlify(im.view( 3, 3,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 8, 3,1,1).tostring()), hex_v);
eq_(hexlify(im.view(13, 3,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 3, 8,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 8, 8,1,1).tostring()), hex_v);
eq_(hexlify(im.view(13, 8,1,1).tostring()), hex_a);
eq_(hexlify(im.view( 3,13,1,1).tostring()), hex_v);
eq_(hexlify(im.view( 8,13,1,1).tostring()), hex_b);
eq_(hexlify(im.view(13,13,1,1).tostring()), hex_v);
def test_rgba_8bui_subquery():
_test_rgba_subquery('rgba_8bui_subquery', '8BUI', 255, 0, 0, 255, 255, 255)
#def test_rgba_16bui_subquery():
# _test_rgba_subquery('rgba_16bui_subquery', '16BUI', 65535, 0, 0, 65535, 65535, 65535)
#def test_rgba_32bui_subquery():
# _test_rgba_subquery('rgba_32bui_subquery', '32BUI')
atexit.register(postgis_takedown)
def enabled(tname):
enabled = len(sys.argv) < 2 or tname in sys.argv
if not enabled:
print "Skipping " + tname + " as not explicitly enabled"
return enabled
if __name__ == "__main__":
setup()
fail = run_all(eval(x) for x in dir() if x.startswith("test_") and enabled(x))
exit(fail)