1

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.

How to calculate the area of polygon overlap in R?

Community
  • 1
  • 1
user3281487
  • 97
  • 1
  • 1
  • 9
  • 1
    Well, `polygon` _"draws the polygons whose vertices are given in x and y"_ and what you need is a `SpatialPolygons` object. Look at the `sp` package a bit more closely. Also, you can use `rgeos::gArea` to calculate the area of the given geometry. It's also a nice idea to include all the necessary `library` calls in a question vs make folks guess at what you did. – hrbrmstr Sep 29 '15 at 01:01
  • +1 to everything @hrbrmstr's comment. For computing the area of polygons in lat-long coordinates, like those you've got, see `areaPolygon()` in the **geosphere** package. – Josh O'Brien Sep 29 '15 at 01:27
  • Thanks so much for the comments, I got it to work! – user3281487 Sep 29 '15 at 16:51
  • than post the code of your solution for the next person and make your question answerd. – and-bri Oct 16 '16 at 04:19

0 Answers0