Calculating area of polygons created from a spatial intersect

Hello all,

I am attempting to automate an analysis that I developed with ArcInfo
using R and the gdal and geos packages (or any other) if possible.

Here is the basic process

I have a shape file (lines) that defines the limits of all of the
projects with each project having a unique identifier.

I have another shape file (polys) that contains total population and
low income population and represent Census block groups.  This shape
file has an area field which has the acreage of the total block group.

Process

Step 1.
I then buffer these project lines to create a second shape file that
represents the 'footprint' of the project. (Creates polys).

Step 2.
In ArcInfo, I perform an intersection of the two shape files
(footprint and census blocks) and this creates a third shape file
which has a unique polygon for every project/census block intersection.

Step 3.
I then perform an area calculation (acres) on this new poly shape file
and use this calculated area divided by the original area of the
associated census block group to apportion the two population datum to
this new polygon.

Step 4.
Finally, I sum the two population datums for each of the projects from
the attribute table of this final shape file.

When I try to replicate the above procedure I run into a problem with
Step 2 when I use what I think is the appropriate command:

gIntersects(buffered_projects, census_blocks, byID=TRUE)

This command is producing a matrix of each project/census block
combination and only providing me a true/false indication.  Is there
any way to replicate the process from ArcInfo that I outlined above
within R?

Walter Anderson