[prev in list] [next in list] [prev in thread] [next in thread] 

List:       r-sig-geo
Subject:    Re: [R-sig-Geo] intersection of polygons
From:       Edzer Pebesma <edzer.pebesma () uni-muenster ! de>
Date:       2013-03-29 14:19:23
Message-ID: 5155A2EB.4070709 () uni-muenster ! de
[Download RAW message or body]



On 03/29/2013 02:20 PM, Ross Ahmed wrote:
> Many thanks Roger. 
> 
> I thought proj4string(spAll) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84") took \
> care of projection, but obviously not. 

No, it only sets the proj4string, which is NA if not set; +proj=longlat
indicates unprojected data, i.e. coordinates reflect longitude and latitude.

> How do I project?

with spTransform in package rgdal - you need to choose a target CRS for
this.

> 
> Ross
> 
> On 29 Mar 2013, at 12:54, Roger Bivand <Roger.Bivand@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@nhh.no
> > 
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma@wwu.de

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


[prev in list] [next in list] [prev in thread] [next in thread]