[R-sig-Geo] [rgeos] A gIntersection operation causes RStudio to crash

[Mathieu Rajerison]

Hi List,


I'm trying to calculate fragstats metrics, in particular fractal dimension
using spatstat package, raster package and rgeos.

Fractal dimension is here applied to buildings, so as to determine
complexity of shapes. I'm calculating it with a moving window (a circular
disc shape).

There's no problem with the 27 first moving windows but my script causes
RStudio to crash when it's about to calculate the intersection of my 28th
disc and my buildings that intersect the latter.

So, I gave you two images of the couple disc + intersecting buildings in
case you'd see the problem,
and a zip containing disc.Rdata + bati.cddtes.Rdaya if yopu'd like to give
it a try

The code:
library(rgeos)

load("bati.cddtes.RData")
load("W.RData")

gIntersection(bati.cddtes, W)


Here are the images:
http://hpics.li/b0a2086 -(zoomed on disc)
http://hpics.li/7e6f780 - (zoomed on buildings)


Thanks for any help!

Mathieu





--
and here is the complete code, for information:

library(spatstat)
library(raster)
library(rgeos)
library(maptools)


# LOAD
load("bati.iris.RData")

proj4string(bati.iris)=CRS("+init=epsg:2154")

# POLYGONS
bati.pol <- bati.iris

# CHECK POLYGON HOLES
bati.pol1 <- slot(bati.pol, "polygons")
bati.pol1 <- lapply(bati.pol1, checkPolygonsHoles)
slot(bati.pol, "polygons") <- bati.pol1
proj4string(bati.pol)=CRS("+init=epsg:2154")

# LINES
bati.line <- as(bati.iris, "SpatialLines")
proj4string(bati.line)=CRS("+init=epsg:2154")


# PARAMETERS
resolution <- 50
radius <- 100

# CREATION DE LA GRILLE RASTER
r <- raster(extent(bati.iris))
res(r) <- resolution
coords <- xyFromCell(r, 1:ncell(r))

## LOOP
#area <- vector(mode="list", length=ncell(r))
#perimeter <- vector(mode="list", length=ncell(r))
seqce <- seq(1:27)
#seqce <- seq(ncell(r))
area <- vector(mode="list", length=length(seqce))
perimeter <- vector(mode="list", length=length(seqce))

#for (i in seq(1:ncell(r))) {
for (i in seq(along=seqce)) {
  print (i)
  # WINDOW
  W <- disc(radius, coords[i,]); W <- as(W, "SpatialPolygons");
proj4string(W)=CRS("+init=epsg:2154")

  # REQUETE SPATIALE
  candidates.pol <- gIntersects(bati.pol, W, byid=TRUE)
  if (any(candidates.pol==TRUE)) {
    # AREA
    res.pol <- gIntersection(bati.pol[candidates.pol[1,],], W)
    area[[i]] <- gArea(res.pol)

    # PERIMETER
    candidates.line <- gIntersects(bati.line, W, byid=TRUE)
    res.line <- gIntersection(bati.line[candidates.line[1,],], W)
    perimeter[[i]] <- gLength(res.line)
  }
}

# CALCULATIONS
# (P9) Fractal Dimension Index
#FRAC[[i]] <- 2*log(0.25*perimeter)/log(area)


  #row.names(PARA[[i]]) <- as.character(i)
  #row.names(FRAC[[i]]) <- as.character(i)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20111017/95c15c44/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fragstats.zip
Type: application/zip
Size: 7647 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20111017/95c15c44/attachment.zip>