Hi List,<br><br><br>I'm trying to calculate fragstats metrics, in particular fractal dimension using spatstat package, raster package and rgeos.<br><br>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).<br>
<br>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.<br><br>So, I gave you two images of the couple disc + intersecting buildings in case you'd see the problem, <br>
and a zip containing disc.Rdata + bati.cddtes.Rdaya if yopu'd like to give it a try <br><br>The code:<br>library(rgeos)<br><br>load("bati.cddtes.RData")<br>load("W.RData")<br><br>gIntersection(bati.cddtes, W)<br>
<br><br>Here are the images:<br><a href="http://hpics.li/b0a2086">http://hpics.li/b0a2086</a> -(zoomed on disc)<br><a href="http://hpics.li/7e6f780">http://hpics.li/7e6f780</a> - (zoomed on buildings)<br><br><br>Thanks for any help!<br>
<br>Mathieu<br><br><br><br><br><br>--<br>and here is the complete code, for information:<br><br>library(spatstat)<br>library(raster)<br>library(rgeos)<br>library(maptools)<br><br><br># LOAD<br>load("bati.iris.RData")<br>
<br>proj4string(bati.iris)=CRS("+init=epsg:2154")<br><br># POLYGONS<br>bati.pol <- bati.iris<br><br># CHECK POLYGON HOLES<br>bati.pol1 <- slot(bati.pol, "polygons")<br>bati.pol1 <- lapply(bati.pol1, checkPolygonsHoles)<br>
slot(bati.pol, "polygons") <- bati.pol1<br>proj4string(bati.pol)=CRS("+init=epsg:2154")<br><br># LINES<br>bati.line <- as(bati.iris, "SpatialLines")<br>proj4string(bati.line)=CRS("+init=epsg:2154")<br>
<br><br># PARAMETERS<br>resolution <- 50<br>radius <- 100<br><br># CREATION DE LA GRILLE RASTER<br>r <- raster(extent(bati.iris))<br>res(r) <- resolution<br>coords <- xyFromCell(r, 1:ncell(r))<br><br>## LOOP<br>
#area <- vector(mode="list", length=ncell(r))<br>#perimeter <- vector(mode="list", length=ncell(r))<br>seqce <- seq(1:27)<br>#seqce <- seq(ncell(r))<br>area <- vector(mode="list", length=length(seqce))<br>
perimeter <- vector(mode="list", length=length(seqce))<br><br>#for (i in seq(1:ncell(r))) {<br>for (i in seq(along=seqce)) {<br>  print (i)<br>  # WINDOW<br>  W <- disc(radius, coords[i,]); W <- as(W, "SpatialPolygons"); proj4string(W)=CRS("+init=epsg:2154")  <br>
  <br>  # REQUETE SPATIALE<br>  candidates.pol <- gIntersects(bati.pol, W, byid=TRUE) <br>  if (any(candidates.pol==TRUE)) {<br>    # AREA<br>    res.pol <- gIntersection(bati.pol[candidates.pol[1,],], W)<br>    area[[i]] <- gArea(res.pol)<br>
    <br>    # PERIMETER<br>    candidates.line <- gIntersects(bati.line, W, byid=TRUE)    <br>    res.line <- gIntersection(bati.line[candidates.line[1,],], W)<br>    perimeter[[i]] <- gLength(res.line)<br>  }<br>
}<br><br># CALCULATIONS<br># (P9) Fractal Dimension Index<br>#FRAC[[i]] <- 2*log(0.25*perimeter)/log(area)<br>  <br>  <br>  #row.names(PARA[[i]]) <- as.character(i)<br>  #row.names(FRAC[[i]]) <- as.character(i)<br>
<br><br>