user21227
user21227

Reputation: 13

Area.poly in R: why Shapefile converts into class "list" instead of "gpc.poly" ?

I am trying to calculate overlap between and two shapefiles using "area.poly" in R. I have use the following coding and my problem is following:

I try to convert the shapefile2 into gpc.poly, the result is that the shapefile gets into a class "list". Therefor my use of area.poly do not work.

Anyone know why the shapefile1 do now want to transform into the gpc.poly class ??

the shapefile contain the following extensions: dbf, prj, qpj, shp, shx

the code I want to use:

-> library(rgeos)

-> p1 <- as(shapefile1, "gpc.poly") -> p2 <- as(shapefile2, "gpc.poly")

-> area.poly(intersect(p1,p2)) get error here!

-> str(shapefile1) result: class(list)

Upvotes: 1

Views: 275

Answers (1)

jlhoward
jlhoward

Reputation: 59375

Why do you want to use intersect(...) and area.poly(...)?? Use gIntersection(...) and gArea(...) instead.

Note that most of this is just to set up the demo. You already have the two shapefiles, apparently, so you would start near the very end.

library(raster)  # for getData(...); just for this example
library(rgdal)   # for spTransform
library(rgeos)   # for gIntersection(...), gArea(...), etc.
# map of the UK - just an example
UK      <- getData("GADM",country="GBR",level=0)
OSGB.36 <- "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs"
UK      <- spTransform(UK,CRS(OSGB.36))   # need projection with meaningful units
# circle with 300km radius from centroid of UK
ctr <- coordinates(gCentroid(UK)) 
th  <- 2*pi*seq(0,.99,by=0.01)
pts <- data.frame(long=ctr[1]+3e5*cos(th),lat=ctr[2]+3e5*sin(th))
pts <- rbind(pts,pts[1,])
circle <-  SpatialPolygons(list(Polygons(list(Polygon(pts)),1)),
                           proj4string=CRS(proj4string(UK)))

## you start here...
plot(UK)
plot(circle,add=T)
sp <- gIntersection(UK,circle)  # calculate intersection
plot(sp,add=T, col="red")
gArea(sp)/1e6                   # calculate area of intersection in sq km
# [1] 149638.8

Upvotes: 0

Related Questions