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

[David Winsemius]

On May 14, 2010, at 10:59 AM, David Winsemius wrote:

>
> 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

And if you use a different function with your original polygon, you  
may get a different, albeit still imperfect, answer:

 > require(sgeostat)
 > my.pol <- data.frame( "x" = c(0,0,100), "y" = c(0,100,0) )
 > in.polygon(
+    x0 = test.bug[,"SILT"],
+    y0 = test.bug[,"CLAY"],
+    x   = my.pol$"x",
+    y   = my.pol$"y"
+ )
[1]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE

-- 
David.

>
> 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
>
> ______________________________________________
> 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