Some points from PointDatasource get lost on reprojection

[artemp]

See attached python script to reproduce. Basically, I create a PointDatasource, add several points with longitude/latitude coordinates to it and attach it to a layer with appropriate projection setting ('+proj=latlon +datum=WGS84'). Then I try to render a map with this layer. Depending on the map projection setting, one or more points get lost.

In an ordinary mercator projection ('+proj=merc +datum=WGS84 +k=1.0 +units=m +over +no_defs'), reproducibly the southernmost point vanishes. In fact, I can work around the problem by setting an invisible 'dummy' point south of the area of interest.

In UTM, additionally the westernmost point disappears.

On the other hand, a map in "Google" spherical mercator projection displays all points. So does a simple latlon map where no reprojection is necessary.

Reading the same points from a shapefile works for all these output projections.

[artemp]

[phispi] I see the same effect using version 6.1 and a GPX file datasource via OGR: As soon as I select an output projection were a reprojection is necessary, much fewer points are there as there should be. There are no problems with lines and reprojection in my case.

I try to setup a debug environment for mapnik so that I can maybe find the problem and supply a patch: Unfortunately, DebuggingMapnik is not very informative ;-)

[artemp]

[phispi] I found the reason :-)

In source:/trunk/src/envelope.cpp@1325#L230 a check is performed whether two envelopes intersect. They should do so, however, for the points that are not shown this function returns false and therefore they are not drawn:

{{{
#!cpp
bool Envelope::intersects(const Envelope &other) const
{
return !(other.minx_>maxx_ || other.maxx_<minx_ ||
other.miny_>maxy_ || other.maxy_<miny_);
}
}}}

The variables for a case where the function returns false where it should not are:

{{{
this:
minx_ 142.4800000000011
miny_ -38.59999999999934
maxx_ 143.30000000001135
maxy_ -38.299999999990874

other:
minx_ 142.48
miny_ -38.38
maxx_ 142.48
maxy_ -38.38
}}}

The term other.maxx_<minx_ gives true but but it should give false because the two values are identical except roundoff errors.

[artemp]

[podolsir] > ... but it should give false because the two values are identical except roundoff errors.
Makes me remember my Intro to Computer Science course: "Never ever compare two floats directly!" ;) Now there's just one problem: how do you tell what is a rounding error and what isn't?

The solution would IMHO be using some kind of tolerance when comparing the values. A fixed tolerance wouldn't probably work very well, though -- a 0.00001 difference doesn't matter in pixmap space or some projected coordinate systems, but does when using geographic coordinates. So it should be relative to the values or something like that (see for example [http://www.mrupp.info/Data/2007floatingcomp.pdf]). Boost has some functions (or function-like classes...) to compare floats safely:

[http://www.boost.org/doc/libs/1_34_1/libs/test/doc/components/test_tools/floating_point_comparison.html]

But it's in the unit test framework which isn't used by mapnik as of now... Maybe there's some other way to integrate it or some kind of such a comparison here.

Unfortunately, I don't really have time now to create a patch and test this at least with the Boost one (and as envelope() is called rather often, I don't know if this would introduce some slowdown to the rendering process, so I would also look at that too). Phispi, maybe you could do that?

[artemp]

[phispi] > Now there's just one problem: how do you tell what is a rounding error and what isn't?

Exactly. The quickest solution would of course be to a comparison that is tolerant to rounding errors. However, I'm not sure that this is the "best" solution in this situation.

One question is where the rounding error comes from. I guess (but I haven't looked it up in the code yet) that the coordinates are projected and projected back - so if that should be true it maybe can be avoided.

An other solution would be to add the tolerance to the bounding box of the point source and maybe there are other possibilities to consider.

I'll assign the ticket to me and play with it - hopefully resulting in my "recommended patch". As I don't have so much time too, it could take a while but when I'm looking at the unclosed tickets for 0.6.2, it doesn't look too urgent for 0.6.2. ;-)

By the way: I think it would make sense to add a test for this kind of problem so that it is detected automatically in the future.

[artemp]

[phispi] Thanks to the help of Gregor, I found a suggestion for a solution sooner as expected :-)

In source:trunk/include/mapnik/feature_style_processor.hpp@1326#L105 the following method was used to determine, whether a point lies within the region that should be drawn or not:

1. The projection of the layer was transformed to the projection of the map
2. The intersection of the two boundaries was calculated
3. The intersection was transformed back to the projection of the layer. The floating point errors introduced in this step caused the problem.

The suggested solution is the following:

1. Take the projection of the map and transform it to the projection of the layer.
2. Determine the intersection here.

This saves one reprojection ''and'' does not introduce the reprojection floating point errors. I upload the short patch. With this patch applied, the pointsymbolizer-test.py script produces valid results where all points are shown. Theoretically this patch should make mapnik even faster ;-)

[artemp]

[phispi] What remains to to is to judge the patch, apply/commit the patch if accepted and to write a test case for this ticket.

[artemp]

[podolsir] Well... While minimizing the number operations like reprojections improves the performance and readability and is therefore worthwhile in its own right, is doesn't solve the underlying core problems that come with floating point numbers -- it only alleviates them in some cases.

Envelopes and coordinates and other floats come to mapnik from all kinds of sources, and they can be subject to any number of operations that can introduce rounding errors and are out of control of mapnik or even the mapnik caller. For example, the envelopes of tiles for the OSM map are computed on the fly from the single bounding box the user inputs. So problems like these will still occur rarely and unpredictably even if they don't anymore in this particular case.

I still think that introducing a proper float comparison is the way to go.

Don't get me wrong: your patch still improves the map processing - but it addresses a different aspect :) As to the patch itself, I looked at it for an half an hour now and checked whether it actually does the same thing as the old code, and I think it's good. One minor thing, though: the variables {{{lx0}}}, {{{ly0}}} and so on seem a little bit misnamed now, as the {{{l}}} stood for layer and you're now storing the ''map'' bbox in them. It is quite irritating when reading the new code, at least for me.

[artemp]

[phispi] > I still think that introducing a proper float comparison is the way to go.

In this case I asked myself: Is it really possible to make a valid error tolerant comparison? The coordinates either are in or out of the envelope. When making the envelope a little big bigger so that points that lie "a little bit" outside the box (due to errors) are still considered to be inside (which would have helped in this ticket case), it might introduce a problem when a point that lies near the border between two envelopes is considered to be within both envelopes that do not intersect.

So I avoided that decision by not introducing calculation errors beforehand. That might not be possible in other cases and some kind or floating point error handling might have to be introduced but for now it could be avoided.

Don't get me wrong: your patch still improves the map processing - but it addresses a different aspect :) As to the patch itself, I looked at it for an half an hour now and checked whether it actually does the same thing as the old code, and I think it's good.

Thank you!

One minor thing, though: the variables {{{lx0}}}, {{{ly0}}} and so on seem a little bit misnamed now, as the {{{l}}} stood for layer and you're now storing the ''map'' bbox in them. It is quite irritating when reading the new code, at least for me.

Ah, OK, sorry, I haven't realized that l stands for layer. In this case, the variables should be renamed of course.

[artemp]

[springmeyer] Hey all, nice discussion and patch! I've been following via rss and just had a chance to look at and test the patch. It does work great for that test case (and we should certainly add this test and many more like it to the source:trunk/tests folder), and it does appear to offer a very small speed bump in rough testing.

I'll try to find time to apply next week.

• dane

[artemp]

[jburgess777] If I remember correctly the original clipping scheme just used a single projection and this had some other issues. The forward + backward scheme was done to fix those.

The way that the data filtering works it is generally OK to return slightly too much data. When fetching lines & polygons we often get back geometries which have parts outside of the requested bbox. We even see objects which are entirely outside the bbox due to the limitations of simple square bbox tests. The rendering code copes with this.

If we multiplied the size of the bbox by 1.00000001 before applying the test then we should get back the lost points. It should not matter that a few points occasionally get returned which really should be outside the box.

[artemp]

[phispi] > If I remember correctly the original clipping scheme just used a single projection and this had some other issues. The forward + backward scheme was done to fix those.

I just looked through the SVN history to see which issues you are referring to. The forward/backward method was introduced in r522 without any comment that this fixes a bug. Before that revision, the layers extend was not clipped to the map's extension at all. Have I overlooked anything?

If we multiplied the size of the bbox by 1.00000001 before applying the test then we should get back the lost points.

I agree - that would solve the problem, too.

Whatever solution is chosen by the mapnik developers, I adapted the pointsymbolizer-test.py script by JRohrer to be included in the test suite. It uses the cairo renderer because here, I can check easily in the SVG string whether the points are present - but maybe a more general solution would be preferable. However, it's much better than nothing. I tested the test: Without patch if fails and with patch it runs through.

The second version of the patch just has different variables names: mx0 instead of lx0 (and so on) as suggested by podolsir.

Thank you all for mapnik by the way - I like it very much!
Philipp

[artemp]

[jburgess777] r522 was done in response to a problem I reported where the extent of a shapefile wraps during projection. See the thread @ https://lists.berlios.de/pipermail/mapnik-devel/2007-September/000243.html

[artemp]

[phispi] > See the thread @ https://lists.berlios.de/pipermail/mapnik-devel/2007-September/000243.html

Thank you for the link! Interesting. Do I understand it correctly that the solution to the problem of 2007 was to clip the layer's extend to the map's bounding box? Therefore the projection/reprojection method was introduced in r522 which solved the problem but introduced "this ticket".

However, with the suggested patch to this ticket, the map is still clipped - so the 2007 problem should still be fixed. But as I didn't completely understand the former problem, I might as well be wrong.

[artemp]

[springmeyer] thanks for the revised patch and test, applied in r1348

[artemp]

[springmeyer] re-opening as it looks like we need to revert this patch.

see: http://www.mail-archive.com/mapnik-devel@lists.berlios.de/msg00547.html

[artemp]

[springmeyer] reverted in r1530 in the 0.7 branch which will be the next release, sorry guys. This fix opens this back up by closes larger problem noted in #486.

see BoundsClipping for more details. In trunk we'll likely apply a new solution to clip the map buffered extent by layer bounds, but for the 0.7 release this will not be fixed.

pushing re-opened ticket to 0.8.0 milestone (current trunk code).

[artemp]

[phispi] Today I played again with this bug because it keeps me from using mapnik. Most of the time I spent trying to compile it (revision r1805) because Debian stable has rather old libraries and mapnik needs new ones ;-) If finally succeeded but running the tests gave me
{{{
Ran 93 tests in 1.562s
FAILED (TODO=6, errors=5, failures=7)
}}}
So I don't know whether the tests fail because of a side-effect caused by mixed libraries on my system or they fail on all systems...

However, this and the #308 ticket seem to be caused by not handling possible errors when using projections like prj_trans.forward(mx0,my0,mz0). The r1348 patch seems to introduce a different behavior when dealing with coordinates that are not possible/valid in one of the used projections. My idea is to modify the reverted r1348 patch to handle coordinates outside the possible coordinate area in a defined way.

To make sure that this does not again make Canada disappear (as mentioned in the conversation cited above) a test case would be useful. I couldn't find one when doing a quick search in the presents tests. Have I overlooked it or do we need to create it?

[artemp]

[springmeyer] Replying to [comment:18 phispi]:

Today I played again with this bug because it keeps me from using mapnik.

Sorry to hear that, we should get on this then. Please also notice the new tickets linked from BoundsClipping: #506, #548, and #549. These all highlight, I think, the need for a patch like the one you've written for this ticket. We just need to have better tests to fully understand the compromises or each approach, as you note below.

Most of the time I spent trying to compile it (revision r1805) because Debian stable has rather old libraries and mapnik needs new ones ;-)

Yes, trunk requires at least Boost 1.41 and icu >= 4.x. See https://trac.mapnik.org/wiki/Mapnik2 for more details.

If finally succeeded but running the tests gave me
{{{
Ran 93 tests in 1.562s
FAILED (TODO=6, errors=5, failures=7)
}}}
So I don't know whether the tests fail because of a side-effect caused by mixed libraries on my system or they fail on all systems...

This is expected with trunk currently. The core library is working fine, but the python api needs updating for the tests to pass (because of changes in the C++ core).

However, this and the #308 ticket seem to be caused by not handling possible errors when using projections like prj_trans.forward(mx0,my0,mz0). The r1348 patch seems to introduce a different behavior when dealing with coordinates that are not possible/valid in one of the used projections. My idea is to modify the reverted r1348 patch to handle coordinates outside the possible coordinate area in a defined way.

One potential approach is to re-apply your patch originally added in r1348, but along with an updated method of constraining the map bounds to valid coordinates. Basically the side affect of r1348, is that moving away from clipping by layer extents meant that one could easily send invalid map extents at zoom level 0/1 (simply by the client application zooming the mapnik map to an area outside of the valid extents of the target projection). #506 starts to get at this, but I think needs to be extended to allow for the user to manually set the max extent (in the case that one of the potentially many map layers has bounds invalid for the target projection).

To make sure that this does not again make Canada disappear (as mentioned in the conversation cited above) a test case would be useful.

Certainly.

I couldn't find one when doing a quick search in the presents tests. Have I overlooked it or do we need to create it?

Yes we need to create it, there is not an existing test that addresses this. Partly because it is a hard thing to capture without visual tests. But either comparing rendered images or using the api to count queried features would be possible (although the latter is not atm possible because changes in the Query class in C++ still need to be exposed in Python (we need to add conversion from python tuples to boost tuples and until then the python query api is broken). You method previously of testing for features rendering to SVG could work as well, but ideally tests can operate without the cairo rendering backend installed.

[artemp]

[springmeyer] I should note that #549 likely has some promising logic, but I've yet to review it.

[artemp]

[olt] I haven't looked at the patches to this ticket, but the cause for missing points is likely the transformation of the bounding boxes. It is not enough to transform only the min and max point. I updated the ticket #594 to take that into account.

[artemp]

[olt] Sorry, I meant #549.

[artemp]

[springmeyer] Thanks everyone for the effort here and especially olt for the nice revised patch over at #549. This ticket, as olt notes, can be solved by the approach at #549 and that is now applied, so closing this ticket.