CODEWITHSUNDEEP
rplotrasterellipse
rar
rar

Reputation: 924

Add an ellipse on raster plot in R

I need to plot an ellipse on a raster plot. I tried plotting raster with simple plot(r1) (r1 is the raster layer) and then add=T for ellipse plotting but it doesn't work. Then I tried axes=F for plotting raster and again tried add=T for ellipse. It still doesn't work.

So I tried converting the ellipse data to dataframe and try adding to raster plot. 

#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)

What I get is; enter image description here

Ideally ps should appear as ellipse but it is a circle. On the other hand if I plot elp (data frame from which ps is created), I get an ellipse.

plot(elp)

enter image description here

Can someone please help with this?

Upvotes: 2

Views: 845

Answers (1)

fdetsch
fdetsch

Reputation: 5308

Be aware that you are not creating an ellipse flattened in horizontal direction at all. As shown in the below console output, the horizontal and vertical extent of ps is almost equal. 

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

I rather assume that the ellipsoid shape of the plot depicted above originates from the sizing of your plotting device. In order to create an ellipse in the first place, you are required to do something like the following (note the difference between the x and y related expansion factors). 

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)))

Once you have a proper ellipse, you may then display it on top of the RasterLayer e.g. by using 

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

spplot

Upvotes: 1

Related Questions