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