I want to calculate the area of overlap between a circle with a given radius on a map and a polygon that partially overlaps with the circle. Here is the map:
map(database="worldHires", xlim = c(-97.1, -80), ylim = c(25, 32.5), col = "grey", fill=TRUE, lwd = 2, bg = "lightblue2", mar=rep(0, 4))
#add Grand Isle
points(-90.0122, 29.2278, pch = 21, col = "blue", bg = "blue", cex=2)
text(-91.5, 28.9, "Grand Isle", cex = 1.2)
Here is the circle with a 100km radius
plotCircle <- function(LonDec, LatDec, Km) {
ER <- 6371 #Mean Earth radius in kilometers. Change this to 3959 and you will have your function working in miles.
AngDeg <- seq(1:360) #angles in degrees
Lat1Rad <- LatDec*(pi/180)#Latitude of the center of the circle in radians
Lon1Rad <- LonDec*(pi/180)#Longitude of the center of the circle in radians
AngRad <- AngDeg*(pi/180)#angles in radians
Lat2Rad <-asin(sin(Lat1Rad)*cos(Km/ER)+cos(Lat1Rad)*sin(Km/ER)*cos(AngRad)) #Latitude of each point of the circle rearding to angle in radians
Lon2Rad <- Lon1Rad+atan2(sin(AngRad)*sin(Km/ER)*cos(Lat1Rad),cos(Km/ER)- sin(Lat1Rad)*sin(Lat2Rad))#Longitude of each point of the circle rearding to angle in radians
Lat2Deg <- Lat2Rad*(180/pi)#Latitude of each point of the circle rearding to angle in degrees (conversion of radians to degrees deg = rad*(180/pi) )
Lon2Deg <- Lon2Rad*(180/pi)#Longitude of each point of the circle rearding to angle in degrees (conversion of radians to degrees deg = rad*(180/pi) )
polygon(Lon2Deg,Lat2Deg,lty=2)
}
plotCircle(-90.0122,29.2278,100)#Plot a dashed circle of 100 km around Grand Isle
This works to plot the circle onto the map. The polygon is here
poly <- readOGR("path/filename")
poly1<- spTransform(poly, CRS("+proj=longlat +datum=WGS84"))
plot(poly1, add=TRUE, col="orange")
This adds the polygon onto the map. I then tried to use
res <-gIntersection(x, poly1)
where x= either plotCircle or a function that combines the map with the circle plotted onto it. I keep getting the following error message:
Error in identical(spgeom1@proj4string, spgeom2@proj4string) : trying to get slot "proj4string" from an object of a basic class ("function") with no slots
I assume that one or the other is not a spatial object, and perhaps it is the circle because I am only drawing it on the map. So I tried to use spTransform
to convert the circle into an SP object but that didn't work. Any ideas on how I could make this happen would be greatly appreciated.
I used the question below as a staring point but it did not solve my problem.