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.
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
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
