Добавьте эллипс на растровый график в R

Question

Мне нужно построить эллипс на растровом графике. Я пробовал рисовать растр с простым plot(r1) (r1 - это растровый слой), а затем add=T для построения эллипса, но он не работает. Затем я попробовал axes=F для печати растра и снова попробовал add=T для эллипса. Это все еще не работает.

Поэтому я попытался преобразовать данные эллипса в dataframe и попробовать добавить к растровому графику.

#Creating a raster 
r <- matrix(sample(1:400),20,20) 
r1<-raster(r) 

#Creating ellipse with given mean and standard deviation values 
theta <- seq(0, 2 * pi, length=(2000)) 
x <- 0.2 - 0.15 * cos(theta) 
y <- 0.5 - 0.15 * sin(theta) 
elp<- cbind(na.omit(x),na.omit(y)) 

#Converting ellipse data frame (elp) to SpatialDataFrame (ps) 
ps <- SpatialPolygons(list(Polygons(list(Polygon(elp)),2))) 

#Plotting raster with ellipse 
plot(r1) 
plot(ps, add=T) 

Что я знаю об этом:

В идеале ps должен выглядеть как эллипс, но это круг. С другой стороны, если я построю elp (кадр данных, из которого создан ps), я получаю эллипс.

plot(elp) 

Может кто-то пожалуйста, помогите с этим?

Answer 1

Помните, что вы не создаете эллипс, сплющенный в горизонтальном направлении. Как показано на приведенном ниже выходе консоли, горизонтальная и вертикальная протяженность ps практически равна.

> diff(range(x)) 
[1] 0.2999998 
> diff(range(y)) 
[1] 0.2999999 

Я скорее предположить, что эллипсоид форма участка, изображенного выше происходит от размеров вашего черчения устройства. Чтобы создать эллипс в первую очередь, вам необходимо сделать что-то вроде следующего (обратите внимание на разницу между связанными факторами расширения x и y).

x <- 0.2 - 0.15 * cos(theta) 
y <- 0.5 - 0.05 * sin(theta) 
elp <- cbind(na.omit(x),na.omit(y)) 
ps <- SpatialPolygons(list(Polygons(list(Polygon(elp)),2))) 

После того, как у вас есть правильный эллипс, вы можете затем отобразить его на верхней части RasterLayer например, с помощью

library(latticeExtra) 
spplot(r1, alpha.regions = 0.5, scales = list(draw = TRUE)) + 
    layer(sp.polygons(ps, col = "black", lwd = 2))