VectorLayer, ExecuteIntersectionQuery, Transformation, is it a bug?

Topics: Algorithms, SharpMap Project
Jan 22, 2015 at 1:02 PM

i had a problem retrieving featuredata from a VectorLayer. i found a solution. it might be considered usefull to integrate the solution in the VectorLayer.cs.

it took me some time to find out where the problem is. i fear it will take some time to show up the problem as well;-)

i had this vectorlayer with exactly 2 points on it. i executed
vectorlayer.ExecuteIntersectionQuery(vectorlayer.envelope, fds)
the resulting table in the featuredataset was empty.

cause of this: my vectorlayer, UTM32 coordinatesystem, works with CoordinateTransformation / ReverseCoordinateTransformation when shown above Openstreetmap.
tracing ExecuteIntersectionQuery the code to transform the box is executed:
#if !DotSpatialProjections
                if (ReverseCoordinateTransformation != null)
                    box = GeometryTransform.TransformBox(box, ReverseCoordinateTransformation.MathTransform);
                    box = GeometryTransform.TransformBox(box, CoordinateTransformation.MathTransform);
                box = GeometryTransform.TransformBox(box, CoordinateTransformation.Target, CoordinateTransformation.Source);
because of the inaccuracy of the transformation, both points were out of box. the inaccuracy is a small one (but enough to cause trouble;-)

box builded of the coordinates of the 2 points:
box = {Env[345198,787433946 : 345657,685866159, 5677448,7220446 : 5680763,68131179]}
box after Transformation:
box = {Env[345198,787433659 : 345657,685865876, 5677448,72204333 : 5680763,68131054]}
my solution was simply extending the box by one meter.
vectorlayer.ExecuteIntersectionQuery(vectorlayer.envelope.ExpandBy(1), fds)
i leave it to the developers of SharpMap to decide wether to implement such expansions after transformations. in my case it was easy to discover the problem. had i done it with a layer of a few thousand points...
Jan 23, 2015 at 9:24 AM
you need to be very careful doing things like this because the value of 1 is relative to the coordinate system you are using. So if you did this using WGS84 coordinate system, you would expand the envelope by 1 degree of latitude/longitude, which would be a big difference.

I think the root of the issue is that points don't really have an envelope because they have no area. I have got round this before by buffering the points a small amount before getting the envelope of the buffered geom.
Jan 27, 2015 at 8:15 AM
you're right, ExpandBy(1) is 1 m using UTM32-coordinates so it worked for me. i didn't take into account that this is different for a coordinatesystem like WGS84. my suggestion to implement some sort of 'correction' wouldn't be trivial. i do not have any idea of possible sideeffects too.

but still it leaves me a bit unsatisfied. if i zoom the map to the layerextend of this 2-point-layer, no point is shown.