[R-sig-Geo] intersection of polygons
Ross Ahmed
rossahmed at googlemail.com
Fri Mar 29 14:20:49 CET 2013
Many thanks Roger.
I thought proj4string(spAll) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84") took care of projection, but obviously not. How do I project?
Ross
On 29 Mar 2013, at 12:54, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> On Fri, 29 Mar 2013, Ross Ahmed wrote:
>
>> Hi all,
>>
>> I have this list of polygons:
>>
>> library(sp)
>> p1 <- Polygon(cbind(c(-1.672487,-1.663689,-1.663185,-1.672258,-1.672487),
>> c(55.58951,55.58948,55.58249,55.58260,55.58951)))
>> p2 <- Polygon(cbind(c(-1.663689, -1.655508, -1.65514, -1.6633, -1.663689),
>> c(55.589485, 55.589554, 55.582396, 55.582428,
>> 55.589485)))
>> p3 <- Polygon(cbind(c(-1.672224, -1.663046, -1.662525, -1.672137,
>> -1.672224),
>> c(55.582142, 55.582044, 55.575139, 55.575348,
>> 55.582142)))
>>
>> sp1 <- Polygons(list(Polygon(p1)),"p1")
>> sp2 <- Polygons(list(Polygon(p2)),"p2")
>> sp3 <- Polygons(list(Polygon(p3)),"p3")
>> spAll <- SpatialPolygons(list(sp1, sp2, sp3))
>> proj4string(spAll) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
>>
>>
>> ?.and also this single polygon:
>>
>> poly1 <- Polygon(cbind(c(-1.666566, -1.659071, -1.658532, -1.666459,
>> -1.666566),
>> c(55.586296, 55.586357, 55.580414, 55.580505,
>> 55.586296)))
>> spoly1 <- Polygons(list(Polygon(poly1)),"poly1")
>> spoly1 <- SpatialPolygons(list(spoly1))
>> proj4string(spoly1) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
>>
>>
>> I need to find the proportion of ?spoly1? that is overlapped by ?sp1?,
>> ?sp2?, ?sp3? and no polygon. The output should be something like:
>>
>> 0.40, 0.30, 0.20, 0.10
>
> library(rgeos)
> gArea(spoly1)
> gArea(gIntersection(spAll, spoly1, byid=TRUE), byid=TRUE)/gArea(spoly1)
>
> but note the warnings. You'll need to project both objects to get a proper measure.
>
> Roger
>
>>
>> Thanks
>> Ross
>>
>>
>>
>> [[alternative HTML version deleted]]
>
> --
> Roger Bivand
> Department of Economics, NHH Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: Roger.Bivand at nhh.no
>
More information about the R-sig-Geo
mailing list