Skip to content

build polygons from ellipses, then sample points within the polygons

5 messages · Fischbach, Anthony, MacQueen, Don, Rolf Turner

#
I have estimated location error ellipses generated from a continuous-time
version of the correlated random walk model for animal telemetry data 
(see http://cran.r-project.org/web/packages/crawl/index.html, functions: 
crwMLE and crwPredict).
Toward this end, I have two questions.

Question 1: Is there a handy function that readily allows creation of
polygons from ellipses?
## toy example
ellipses<-data.frame(x=c(251753, 252966, 254179), 
						y=c(343394, 343315, 343237), 
						majorAxis=c(0.0, 2753.5, 2114.8), 
						minorAxis=c(0.0,  2754.1, 2118.0),
						angle=c(0,0,0)) ## dataFrame with Ellipse information
						## in this case the majorAxis is assumed to lie along the x axis
                                                ## x, y define the centroid
of the ellipse.
coordinates(ellipses) <- x~y  ## cast to a spatialPointsDataFrame
ellipses.poly <-
Wonder.function.that.Makes.spatialPolygons.from.ellipses(ellipses)

Question 2:  I know that there is a handy means to randomly generate
locations from within a SpatialPolygon.  Please help me find it.



-----
Tony Fischbach, Wildlife Biologist
Walrus Research Program
Alaska Science Center
U.S. Geological Survey
4210 University Drive
Anchorage, AK 99508-4650

AFischbach at usgs.gov
http://alaska.usgs.gov/science/biology/walrus
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/build-polygons-from-ellipses-then-sample-points-within-the-polygons-tp7585922.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
1 day later
#
If you have the x,y coordinates of an ellipse, you can create a
SpatialPolygons or SpatialPolygonsDataFrame following the example in ?over
(package sp). 

For sampling, see ?spsample

-Don
#
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?



-----
Tony Fischbach, Wildlife Biologist
Walrus Research Program
Alaska Science Center
U.S. Geological Survey
4210 University Drive
Anchorage, AK 99508-4650

AFischbach at usgs.gov
http://alaska.usgs.gov/science/biology/walrus
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/build-polygons-from-ellipses-then-sample-points-within-the-polygons-tp7585922p7585940.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
#
On 13/03/14 06:55, Anthony Fischbach wrote:
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
6 days later
#
Thank you Rolf.
Below is how I built it in the end.
I had to take measures to ensure that the ellipse polygon closed back on
itself.
Roger:  Would it be appropriate to include this function in sp or an
associated contributed package.
It seems that building spatialPolygons from ellipse specifications would be
a task required of many sp users.

## modified code from Rolf Turner
require(sp)
require(spatstat)
ellipseSP <- function (x0, y0, a, b, phi=pi/2, npts = 100, ...) {? ##
Returns SpatialPolygons ellipse
											      ## phi defaults to 90 degrees 
                                                                                             
## whereby a lies in the N-S and
                                                                                             
## b in the E-W directions
	theta <- seq(0, 2 * pi, length = npts)??????
	xh <- ?a * cos(theta)? ?????
	yh <- ?b * sin(theta)?  ?????
	co <- cos(phi)? ?????
	si <- sin(phi)? ? ???
	x ?<- x0 + co*xh - si*yh
	x[npts]<-x[1]	## force closure of polygon
	y ?<- x0 + si*xh + co*yh??
	y[npts]<-y[1]	## force closure of polygon
	return(SpatialPolygons(list(Polygons(list(with(list(x = x, y =
y),Polygon(cbind(x,y)))),1)))?)
}?  
## Example plot
eSP <- ellipseSP(x0=0, y0=0, a=8, b=4)? 
plot(eSP)?
##use sp::spsample to randomly select points from within the ellipse 
points(spsample(eSP, n=100, type="random"), col='red', pch=19, cex=.25)  
#



-----
Tony Fischbach, Wildlife Biologist
Walrus Research Program
Alaska Science Center
U.S. Geological Survey
4210 University Drive
Anchorage, AK 99508-4650

AFischbach at usgs.gov
http://alaska.usgs.gov/science/biology/walrus
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/build-polygons-from-ellipses-then-sample-points-within-the-polygons-tp7585922p7585978.html
Sent from the R-sig-geo mailing list archive at Nabble.com.