[R-sig-Geo] Calculate Area of Intersection

[St John Brown]

Hi All,

I would like calculate the area of the intersection between a given city in the US and a 'circle' of a certain radius and center location. Below is my implementation.

I think I have done this right. I would like to double check to see if I have made some glaring error, or if there was a better way to implement the task.

I really appreciate any help you may provide.

Thanks!

library(sp)
library(rgdal)
library(maps)
library(rgeos)
library(mapproj)
library(geosphere)
library(maptools)
library(UScensus2010)
library(UScensus2010cdp)

make_circle_region = function(circle_center, circle_radius, earth_radius=3958.76) { #radius in miles
  circle_region = data.frame(matrix(ncol = 2, nrow = 361))
  colnames(circle_region) = c('longitude', 'latitude')
  for (bearing in 0:360){
    circle_region[bearing+1,] = destPoint(p=circle_center,
b=bearing, d=circle_radius, r=earth_radius)
  }
  return(circle_region)
}

map("state", "texas")

austin = city("austin", "texas")
plot(austin, col="red", add=TRUE)

#Plot circle region around an zip code 78721 with radius = 10 miles
circle_center = c(-97.68665, 30.27293 ) #lon & lat of zip code 78721
circle_radius = 10 #units in miles (same units as earth_radius units)
truck_circle = make_circle_region(circle_center, circle_radius)
polygon(mapproject(truck_circle$longitude, truck_circle$latitude),col=rgb(0,1,0,0.5))

#Make circle a spatial object
truck_sp =
SpatialPolygons(list(Polygons(list(Polygon(truck_circle, hole=FALSE)), 1)))
proj4string(truck_sp) = CRS("+proj=longlat +datum=WGS84") #explicitly define the input projections using PROJ4 syntax

# Project to Albers Equal Area, with parameters set for North America
proj.string.alb = CRS("+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=WGS84  +units=km +no_defs")
austin_alb = spTransform(austin, CRS=proj.string.alb)
truck_alb = spTransform(truck_sp, CRS=proj.string.alb)

#Calculate area of intersection
myInt = gIntersection(austin_alb, truck_alb)
sprintf("The intersection is %.2f square km", gArea(myInt))