[R-sig-Geo] build polygons from ellipses, then sample points within the polygons

[Rolf Turner]

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