Given an SpatialPolygons DF with i.e. 2 intersecting polygons, A and
,B, is there a function that would split the original polygons into "onlyA", "onlyB" and "AintersectingB" polygons? Thanks Agus _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
On 29/06/18 00:51, Agustin Lobo wrote: > Given an SpatialPolygons DF with i.e. 2 intersecting polygons, A and > ,B, is there a function > that would split the original polygons into "onlyA", "onlyB" and > "AintersectingB" polygons? You can do this easily in spatstat by converting your polygons to "owin" objects. cheers, Rolf Turner -- Technical Editor ANZJS Department of Statistics University of Auckland Phone: +64-9-373-7599 ext. 88276 _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Mimicking your vignette:
require(rgdal) require(raster) require(spatstat) require(rgeos) require(maptools) dzfile <-"https://www.dropbox.com/s/xxujcwqy3ec21sa/TDM1_DEM__04_S63W062_epsg32720.zip?dl=0" download.file(dzfile,"TDM1_DEM__04_S63W062_epsg32720.zip",method="wget") unzip("TDM1_DEM__04_S63W062_epsg32720.zip") v <- readOGR(dsn="TDM1_DEM__04_S63W062_epsg32720", layer="TDM1_DEM__04_S63W062_epsg32720", stringsAsFactors = FALSE) plot(v) p <- slot(v, "polygons") p2 <- lapply(p, function(x) { SpatialPolygons(list(x)) }) w <- lapply(p2, as.owin) #maptools required te <- tess(tiles=w) class(te) plot.tess(te,do.labels=TRUE) But this is still the original polygons, not the mosaic of polygon parts I'm looking for. Would I need to go doing all possible intersections and then adding the "A not B" and "B not A" parts for all possible pairs? Or is there a function in spatstat to convert "te" into a tesselation of adjacent polygons? Thanks On Thu, Jun 28, 2018 at 2:02 PM, Rolf Turner <[hidden email]> wrote: > > On 29/06/18 00:51, Agustin Lobo wrote: > >> Given an SpatialPolygons DF with i.e. 2 intersecting polygons, A and >> ,B, is there a function >> that would split the original polygons into "onlyA", "onlyB" and >> "AintersectingB" polygons? > > > You can do this easily in spatstat by converting your polygons to "owin" > objects. > > cheers, > > Rolf Turner > > -- > Technical Editor ANZJS > Department of Statistics > University of Auckland > Phone: +64-9-373-7599 ext. 88276 > > _______________________________________________ > R-sig-Geo mailing list > [hidden email] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Agustin: I have attached a function that I think does what you want. It returns either a tessellation (if you leave tess=TRUE) or a list of owin objects otherwise. I *think* that attachments with .R extensions make it through to the list. The attachment should at least get through to Agustin (who is the person of central interest!). It produces all non-empty intersections between sets of tiles in the tessellation argument "ttt". If singles=TRUE (the default) these "intersections" include the tiles themselves. Otherwise at least two tiles are involved in each intersection. E.g. x1 <- allIntersections(te) # a tessellation x2 <- allIntersections(te,tess=FALSE) # a list of windows x3 <- allIntersections(te,singles=FALSE) # tiles themselves omitted. Note that the code includes a work-around for an infelicity in in intersect.owin(). This infelicity will very likely be fixed in a future release of spatstat. cheers, Rolf -- Technical Editor ANZJS Department of Statistics University of Auckland Phone: +64-9-373-7599 ext. 88276 On 02/07/18 23:35, Agustin Lobo wrote: > Mimicking your vignette: > > require(rgdal) > require(raster) > require(spatstat) > require(rgeos) > require(maptools) > > dzfile <-"https://www.dropbox.com/s/xxujcwqy3ec21sa/TDM1_DEM__04_S63W062_epsg32720.zip?dl=0" > download.file(dzfile,"TDM1_DEM__04_S63W062_epsg32720.zip",method="wget") > unzip("TDM1_DEM__04_S63W062_epsg32720.zip") > v <- readOGR(dsn="TDM1_DEM__04_S63W062_epsg32720", > layer="TDM1_DEM__04_S63W062_epsg32720", stringsAsFactors = FALSE) > plot(v) > p <- slot(v, "polygons") > p2 <- lapply(p, function(x) { SpatialPolygons(list(x)) }) > w <- lapply(p2, as.owin) #maptools required > > te <- tess(tiles=w) > class(te) > plot.tess(te,do.labels=TRUE) > > But this is still the original polygons, > not the mosaic of polygon parts I'm looking for. > Would I need to go doing all possible intersections and then > adding the "A not B" and "B not A" parts for all possible pairs? > > Or is there a function in spatstat to convert "te" into a tesselation > of adjacent polygons? > > Thanks > > > > On Thu, Jun 28, 2018 at 2:02 PM, Rolf Turner <[hidden email]> wrote: >> >> On 29/06/18 00:51, Agustin Lobo wrote: >> >>> Given an SpatialPolygons DF with i.e. 2 intersecting polygons, A and >>> ,B, is there a function >>> that would split the original polygons into "onlyA", "onlyB" and >>> "AintersectingB" polygons? >> >> >> You can do this easily in spatstat by converting your polygons to "owin" >> objects. >> >> cheers, >> >> Rolf Turner >> >> -- >> Technical Editor ANZJS >> Department of Statistics >> University of Auckland >> Phone: +64-9-373-7599 ext. 88276 >> >> _______________________________________________ >> R-sig-Geo mailing list >> [hidden email] >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo allIntersections.R (1K) Download Attachment |
Agustin Lobo writes:
>>> Given an SpatialPolygons DF with i.e. 2 intersecting polygons, A and B, >>> is there a function that would split the original polygons into "onlyA", "onlyB" and >>> "AintersectingB" polygons? First convert the SpatialPolygons to spatstat 'owin' objects using as.owin. If there are two polygons A and B, you can just use intersect.owin and setminus.owin to get what you want: AandB <- intersect.owin(A, B ) AnotB <- setminus.owin(A, B ) BnotA <- setminus.owin(B, A) If there are several windows, make a list P containing the windows. Then call Z <- kaleido(P) plot(Z, do.col=TRUE) using the following code. The result Z is a tessellation. Extract the individual pieces as tiles(Z). kaleido <- function(P) { U <- union.owin(as.solist(P)) V <- lapply(P, onesplit, U=U) Z <- Reduce(intersect.tess, V) return(Z) } onesplit <- function(X, U) tess(tiles=list(A=X, NotA=setminus.owin(U, X)), W=U) Prof Adrian Baddeley DSc FAA John Curtin Distinguished Professor Department of Mathematics and Statistics Curtin University, Perth, Western Australia [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Free forum by Nabble | Edit this page |