0

Aim

Calculating the area of overlap of two standard deviational ellipses in R and work out their intersection and union.

Tool to be used

calc_sde() from the 'aspace' package: calc_sde()

Question

Given some of the values returned by the function (see below), is it possible to calculate the area of the ellipses and the extent of their overlap? In particular, I would like to be able to work out the intersection and union between the two polygons. Not an expert in trigonometry and maths, so I am wondering which of the following values returned by the aboventioned function can be put to work and how:

CENTRE.x    
X-coordinate of the centre

CENTRE.y    
Y-coordinate of the centre

Sigma.x 
Half-length of axis along x-axis

Sigma.y 
Half-length of axis along y-axis


Theta   
Rotation angle in degrees

Eccentricity    
A measure of eccentricity (i.e., the flatness of the ellipse)

Area.sde    
Area of the SDE

TanTheta    
Trigonometric result

SinTheta    
Trigonometric result

CosTheta    
Trigonometric result

SinThetaCosTheta    
Trigonometric result

Sin2Theta   
Trigonometric result

Cos2Theta   
Trigonometric result

ThetaCorr   
Corrected theta angle for rotation of major axis from north
2

Lets make two sample objects, from the first and second five points in the sample data in the example for calc_sde:

z1 = calc_sde(id=1,calccentre=TRUE,weighted=FALSE, points=activities[1:5,], verbose=TRUE)
z2 = calc_sde(id=1,calccentre=TRUE,weighted=FALSE, points=activities[6:10,], verbose=TRUE)

That function is pretty appaling and I'd suggest you don't use it - it writes a file, it assigns objects in the global environment, and it prints that list out. All things that R functions really shouldn't do. Never mind. It returns a data frame that is the points of the ellipses, so we have the points in the ellipses now.

Using sf package functions we can create polygons:

library(sf)
s1 = st_polygon(list(as.matrix(z1[,2:3])))
s2 = st_polygon(list(as.matrix(z2[,2:3])))

And we can intersect them:

s12intersect = st_intersection(s1,s2)

plot(s1)
plot(s2, add=TRUE)
plot(s12intersect, add=TRUE, col="red")

enter image description here

And print the areas if we want that as well:

st_area(s1)
st_area(s2)
st_area(s12intersect)
2
  • thank you for your suggestion, which worked indeed. Yes, I agree with you about that package. Do you think that the same could be accomplished without using it? Thanks for any further insight.
    – NewAtGis
    May 30 '20 at 16:27
  • I won't judge the whole package on that one function. First, its possible cleaner functions exist in the package - I've not looked. Secondly I'd suggest the author how this could be cleaner. Thirdly, if the author is responsive, I'd hack the function into a cleaner version and some support functions. I don't know how old or supported this package is or if similar functionality is available elsewhere.
    – Spacedman
    May 30 '20 at 17:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.