[R] point.in.polygon() in sp package: accuracy problems?
David Winsemius
dwinsemius at comcast.net
Fri May 14 17:38:37 CEST 2010
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
More information about the R-help
mailing list