improve the map.query_point/query_map_point implementation, now throwing for invalid coords or projection transformations

This commit is contained in:
Dane Springmeyer 2012-04-13 11:28:30 -07:00
parent a4b3daf181
commit 5767c65470
4 changed files with 160 additions and 67 deletions

View file

@ -9,6 +9,8 @@ For a complete change history, see the SVN log.
## Mapnik 2.1.0
- Improved error feedback for invalid values passed to map.query_point
- Fixed rendering of thin svg lines (#1129)
- Improved logging/debugging system with release logs and file redirection (#937 and partially #986, #467)

View file

@ -564,85 +564,66 @@ CoordTransform Map::view_transform() const
featureset_ptr Map::query_point(unsigned index, double x, double y) const
{
if ( index< layers_.size())
if (!current_extent_.valid())
{
throw std::runtime_error("query_point: map extent is not intialized, you need to set a valid extent before querying");
}
if (!current_extent_.intersects(x,y))
{
throw std::runtime_error("query_point: x,y coords do not intersect map extent");
}
if (index < layers_.size())
{
mapnik::layer const& layer = layers_[index];
try
mapnik::datasource_ptr ds = layer.datasource();
if (ds)
{
double z = 0;
mapnik::projection dest(srs_);
mapnik::projection source(layer.srs());
proj_transform prj_trans(source,dest);
prj_trans.backward(x,y,z);
double minx = current_extent_.minx();
double miny = current_extent_.miny();
double maxx = current_extent_.maxx();
double maxy = current_extent_.maxy();
prj_trans.backward(minx,miny,z);
prj_trans.backward(maxx,maxy,z);
double tol = (maxx - minx) / width_ * 3;
mapnik::datasource_ptr ds = layer.datasource();
if (ds)
double z = 0;
if (!prj_trans.equal() && !prj_trans.backward(x,y,z))
{
throw std::runtime_error("query_point: could not project x,y into layer srs");
}
// TODO - pass tolerance to features_at_point as well
featureset_ptr fs = ds->features_at_point(mapnik::coord2d(x,y));
if (fs)
{
mapnik::box2d<double> map_ex = current_extent_;
if (maximum_extent_) {
map_ex.clip(*maximum_extent_);
}
if (!prj_trans.backward(map_ex,PROJ_ENVELOPE_POINTS))
{
std::ostringstream s;
s << "query_point: could not project map extent '" << map_ex
<< "' into layer srs for tolerance calculation";
throw std::runtime_error(s.str());
}
double tol = (map_ex.maxx() - map_ex.minx()) / width_ * 3;
MAPNIK_LOG_DEBUG(map) << "map: Query at point tol=" << tol << "(" << x << "," << y << ")";
featureset_ptr fs = ds->features_at_point(mapnik::coord2d(x,y));
if (fs)
return boost::make_shared<filter_featureset<hit_test_filter> >(fs,
hit_test_filter(x,y,tol));
return boost::make_shared<filter_featureset<hit_test_filter> >(fs,
hit_test_filter(x,y,tol));
}
}
catch (...)
{
MAPNIK_LOG_ERROR(map) << "Exception caught in \"query_point\"";
}
}
else
{
std::ostringstream s;
s << "Invalid layer index passed to query_point: '" << index << "'";
if (layers_.size() > 0) s << " for map with " << layers_.size() << " layers(s)";
else s << " (map has no layers)";
throw std::out_of_range(s.str());
}
return featureset_ptr();
}
featureset_ptr Map::query_map_point(unsigned index, double x, double y) const
{
if ( index< layers_.size())
{
mapnik::layer const& layer = layers_[index];
CoordTransform tr = view_transform();
tr.backward(&x,&y);
try
{
mapnik::projection dest(srs_);
mapnik::projection source(layer.srs());
proj_transform prj_trans(source,dest);
double z = 0;
prj_trans.backward(x,y,z);
double minx = current_extent_.minx();
double miny = current_extent_.miny();
double maxx = current_extent_.maxx();
double maxy = current_extent_.maxy();
prj_trans.backward(minx,miny,z);
prj_trans.backward(maxx,maxy,z);
double tol = (maxx - minx) / width_ * 3;
mapnik::datasource_ptr ds = layer.datasource();
if (ds)
{
MAPNIK_LOG_DEBUG(map) << "map: Query at point tol=" << tol << "(" << x << "," << y << ")";
featureset_ptr fs = ds->features_at_point(mapnik::coord2d(x,y));
if (fs)
return boost::make_shared<filter_featureset<hit_test_filter> >(fs,
hit_test_filter(x,y,tol));
}
}
catch (...)
{
MAPNIK_LOG_ERROR(map) << "Exception caught in \"query_map_point\"";
}
}
return featureset_ptr();
CoordTransform tr = view_transform();
tr.backward(&x,&y);
return query_point(index,x,y);
}
Map::~Map() {}

View file

@ -0,0 +1,107 @@
#!/usr/bin/env python
from nose.tools import *
from utilities import execution_path
from copy import deepcopy
import os, mapnik
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('.'))
# map has no layers
@raises(IndexError)
def test_map_query_throw1():
m = mapnik.Map(256,256)
m.zoom_to_box(mapnik.Box2d(-1,-1,0,0))
m.query_point(0,0,0)
# only positive indexes
@raises(IndexError)
def test_map_query_throw2():
m = mapnik.Map(256,256)
m.query_point(-1,0,0)
# map has never been zoomed (nodata)
@raises(RuntimeError)
def test_map_query_throw3():
m = mapnik.Map(256,256)
m.query_point(0,0,0)
# map has never been zoomed (even with data)
@raises(RuntimeError)
def test_map_query_throw4():
m = mapnik.Map(256,256)
mapnik.load_map(m,'../data/good_maps/agg_poly_gamma_map.xml')
m.query_point(0,0,0)
# invalid coords in general (do not intersect)
@raises(RuntimeError)
def test_map_query_throw5():
m = mapnik.Map(256,256)
mapnik.load_map(m,'../data/good_maps/agg_poly_gamma_map.xml')
m.zoom_all()
m.query_point(0,9999999999999999,9999999999999999)
# invalid coords for back projecting
@raises(RuntimeError)
def test_map_query_throw6():
m = mapnik.Map(256,256)
mapnik.load_map(m,'../data/good_maps/merc2wgs84_reprojection.xml')
wgs84_bounds = mapnik.Box2d(-180,-90,180,90)
m.maximum_extent = wgs84_bounds
m.zoom_all()
m.query_point(0,-180,-90)
def test_map_query_works1():
m = mapnik.Map(256,256)
mapnik.load_map(m,'../data/good_maps/wgs842merc_reprojection.xml')
merc_bounds = mapnik.Box2d(-20037508.34,-20037508.34,20037508.34,20037508.34)
m.maximum_extent = merc_bounds
m.zoom_all()
fs = m.query_point(0,-11012435.5376, 4599674.6134) # somewhere in kansas
feat = fs.next()
eq_(feat.attributes['NAME_FORMA'],u'United States of America')
def test_map_query_works2():
m = mapnik.Map(256,256)
mapnik.load_map(m,'../data/good_maps/merc2wgs84_reprojection.xml')
wgs84_bounds = mapnik.Box2d(-179.999999975,-85.0511287776,179.999999975,85.0511287776)
m.maximum_extent = wgs84_bounds
# caution - will go square due to evil aspect_fix_mode backhandedness
m.zoom_all()
#mapnik.render_to_file(m,'works2.png')
# valid that aspec_fix_mode modified the bbox
eq_(m.envelope(),mapnik.Box2d(-179.999999975,-179.999999975,179.999999975,179.999999975))
fs = m.query_point(0,-98.9264, 38.1432) # somewhere in kansas
feat = fs.next()
eq_(feat.attributes['NAME'],u'United States')
def test_map_query_in_pixels_works1():
m = mapnik.Map(256,256)
mapnik.load_map(m,'../data/good_maps/wgs842merc_reprojection.xml')
merc_bounds = mapnik.Box2d(-20037508.34,-20037508.34,20037508.34,20037508.34)
m.maximum_extent = merc_bounds
m.zoom_all()
fs = m.query_map_point(0,55,100) # somewhere in middle of us
feat = fs.next()
eq_(feat.attributes['NAME_FORMA'],u'United States of America')
def test_map_query_in_pixels_works2():
m = mapnik.Map(256,256)
mapnik.load_map(m,'../data/good_maps/merc2wgs84_reprojection.xml')
wgs84_bounds = mapnik.Box2d(-179.999999975,-85.0511287776,179.999999975,85.0511287776)
m.maximum_extent = wgs84_bounds
# caution - will go square due to evil aspect_fix_mode backhandedness
m.zoom_all()
# valid that aspec_fix_mode modified the bbox
eq_(m.envelope(),mapnik.Box2d(-179.999999975,-179.999999975,179.999999975,179.999999975))
fs = m.query_map_point(0,55,100) # somewhere in middle of us
feat = fs.next()
eq_(feat.attributes['NAME'],u'United States')
if __name__ == "__main__":
setup()
[eval(run)() for run in dir() if 'test_' in run]

View file

@ -69,17 +69,20 @@ def test_dataraster_query_point():
_map = mapnik.Map(256,256, srs)
_map.layers.append(lyr)
# point inside raster extent with valid data
x, y = 427417, 4477517
x, y = 556113.0,4381428.0 # center of extent of raster
_map.zoom_all()
features = _map.query_point(0,x,y).features
assert len(features) == 1
feat = features[0]
center = feat.envelope().center()
assert center.x==x and center.y==y, center
value = feat['value']
assert value == 21.0, value
assert value == 18.0, value
# point outside raster extent
# point inside map extent but outside raster extent
current_box = _map.envelope()
current_box.expand_to_include(-427417,4477517)
_map.zoom_to_box(current_box)
features = _map.query_point(0,-427417,4477517).features
assert len(features) == 0