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

