r - 计算特定空间区域中置信椭圆的面积

Question

我想知道是否有人知道如何计算我的置信椭圆内的蓝色阴影区域。任何建议或值得一看的地方都将不胜感激。另外，我希望找到一个通用公式，因为椭圆在应用中不一定必须位于该区域(即椭圆可能更大)。如果有帮助的话，这是我的图片代码:

library(car)

x = c(7,4,1)
y = c(4,6,7)

plot(x,y,xlim=c(0,10),ylim=c(0,10))
rect(x,y,max(x)+100000,max(y)+100000,col="lightblue",border=NA)
points(x,y,col="red",pch=19)

ellipse(center=c(3.5,5),shape=matrix(c(1,.5,.5,2),nrow=2),radius=3,col="green")

Answer 1

如果您可以将椭圆和蓝色区域转换为 SpatialPolygons 对象，那么使用 rgeos 包中的函数来计算它们的面积(面积)就轻而易举了它们的交集，以及所有其他类型的有趣的量。

不幸的是，以正确的形式获取对象需要一些繁重的工作:

library(car)
library(sp)
library(rgeos)

## Function for creating a SpatialPolygons object from data.frame of coords
xy2SP <- function(xy, ID=NULL) {
    if(is.null(ID)) ID <- sample(1e12, size=1)
    SpatialPolygons(list(Polygons(list(Polygon(xy)), ID=ID)),
                    proj4string=CRS("+proj=merc"))
}

## Ellipse coordinates
plot.new()  # Needed by ellipse() 
ell <- ellipse(center=c(3.5,5),shape=matrix(c(1,.5,.5,2),nrow=2),radius=3)
dev.off()   # Cleaning up after plot.new()

## Three rectangles' coordinates in a length-3 list
x <- c(7,4,1)
y <- c(4,6,7)
mx <- max(x) + 1e6
my <- max(y) + 1e6
rr <- lapply(1:3, function(i) {
    data.frame(x = c(x[i], x[i], mx, mx, x[i]),
               y = c(y[i], my, my, y[i], y[i]))
})


## Make two SpatialPolygons objects from ellipse and merged rectangles
ell <- xy2SP(ell)
rrr <- gUnionCascaded(do.call(rbind, lapply(rr, xy2SP)))

## Find area of their intersection
gArea(gIntersection(ell, rrr))
# [1] 10.36296