On 13/03/14 06:55, Anthony Fischbach wrote:
What I am looking for is a means to generate spatialPolygons from the ellipse
specifications (x,y or central coordinate, axis-a and axis-b length, along
with implied orientation of the two axes, which in my case are N-S, E-W).
Is there a handy function enabling this?
Easy enough to write one:
ellipse <- function (x0, y0, a, b, phi, npts = 200, plotit = FALSE,
add = FALSE, phil = FALSE, ...) {
# Build "horizontal" ellipse centred at 0:
theta <- seq(0, 2 * pi, length = npts)
xh <- a * cos(theta)
yh <- b * sin(theta)
# Rotate through angle phi and shift centre:
co <- cos(phi)
si <- sin(phi)
x <- x0 + co*xh - si*yh
y <- x0 + si*xh + co*yh
# Plot the result if required:
if (plotit) {
if (add) {
if (phil)
polygon(x, y, ...)
else lines(x, y, ...)
}
else {
if (phil) {
plot(x, y, type = "n", ...)
polygon(x, y, ...)
}
else plot(x, y, type = "l", ...)
}
return(invisible(list(x = x, y = y)))
}
list(x = x, y = y)
}
This creates an ellipse in the form of a *list*. If you a
SpatialPolgons object you will have to convert it. E.g.:
A <- ellipse(3,7,5,2,pi/4)
B <- with(A,Polygon(cbind(x,y)))
C <- Polygons(list(B),1)
D <- SpatialPolygons(list(C))
plot(D)
The foregoing conversion is a result of trial-and-error mucking around
on my part. It can probably be done in a much sexier manner. I don't
really know what I'm doing with these incomprehensible S4 classes and
methods. Others may advise you as to how to do it properly.
Or you could use the spatstat package:
EW <- owin(poly=A)
X <- runifpoint(50,win=EW) # 50 points distributed uniformly
# random on the interior of the ellipse.
plot(X)
cheers,
Rolf Turner
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
