[R] point.in.polygon() in sp package: accuracy problems?

[David Winsemius]

On May 14, 2010, at 4:50 AM, Julien Moeys wrote:

> Dear list:
>
> I encountered some problems using the function point.in.polygon() of  
> the sp package, when trying to determine whether some points lye  
> inside, outside, on the border or on a vertice of a polygon.

Any time that you are asking questions about what should happen in a  
theoretically perfect, mathematically correct universe, you need to  
carefully what will happen when you map that question over to a  
computationally discrete and approximate computer reality. R is not a  
symbolic algebra system.

>
> I have a list of point I know should lye right on the border of a  
> polygon, but some of them are not classified as such by  
> point.in.polygon() (see the example code below).
>
> To make a long story short I am working on ternary variables that  
> sum to 100% (soil texture data), that can be plotted on a ternary  
> diagram (soil texture diagram) & classified according to  
> subdivisions of the diagram (soil texture class). For a given data  
> point, when one or more of the the 3 variable is equal to 0, and  
> when their sum is 100%, I know the point lies on the border of the  
> ternary diagram. Here I got a case where some of these points are  
> NOT classified as "on the border" of the polygon that defines the  
> limits of my diagram (by point.in.polygon).
>
> The example below shows a simple example with one polygon with a  
> right angle (easy to plot).
>
> I suspect this is a problem of accuracy (as I already got similar  
> problems with a code I made, that proved to be due to accuracy  
> problems arising from from trigonometric calculations), but I would  
> like to know if I am right or if this is another problem (maybe in  
> my own code)... I must say that with 'real & uncertain soil data', I  
> don't really care to know if a point is right on the border of a  
> polygon, but it is good to know what cause the problem!
>
> Thanks for any help
>
> Julien
>
> # --- --- --- ---
> require(sp)
>
> # Dummy ternary data to be plotted & classified as in or out a polygon
> # All should be on the border
> test.bug <- data.frame(
>    #              1      2      3      4      5      6      7      8
>    "CLAY" = c(10.50, 11.05, 20.91, 20.95, 20.00, 00.00, 01.96, 01.24),
>    "SILT" = c(00.00, 00.00, 79.09, 79.05, 80.00, 02.00, 00.00, 00.00),
>    "SAND" = c(89.50, 88.95, 00.00, 00.00, 00.00, 98.00, 98.04, 98.76)
> )   #
>
> # ACoordinates of the polygon
> my.pol <- data.frame( "x" = c(0,0,100), "y" = c(0,100,0) )
>
> # Set up plot scene
> plot( x = my.pol$"x", y = my.pol$"y" )
>
> # Draw the polygon
> polygon(
>    x      = my.pol$"x",
>    y      = my.pol$"y",
>    border = "red",
>    col    = "lightgray"
> )   #
>
> # Draw the point
> points(
>    x   = test.bug[,"SILT"],
>    y   = test.bug[,"CLAY"],
>    pch = 2
> )   #
>
> # Classify the points
> point.in.polygon(
>    point.x = test.bug[,"SILT"],
>    point.y = test.bug[,"CLAY"],
>    pol.x   = my.pol$"x",
>    pol.y   = my.pol$"y"
> )   #
> # [1] 2 2 0 1 2 2 2 2
> # All the points should be 2

And if you widen the polygon by 0.01 at each vertex you get most of  
the points seen as inside and one (the one that was earlier  
outside( now seen on the border:

my.pol <- data.frame( "x" = c(0-0.01, 0-0.01, 100+0.01), "y" =  
c(0-0.01, 100+0.01, 0-0.01) )
(rest of code the same)
[1] 1 1 2 1 1 1 1 1

Now:
Read FAQ 7.31

-- 
David.
>
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius, MD
West Hartford, CT