Hi,

Thanks for the help of Roger Bivand that previously gave me a method on how
to conveniently compute intersections between objets using rgeos and the
help of Colin Rundel, all works fine.

- For those who are interested, here is the result of fractal dimension
index calculation:
http://hpics.li/2e2ff78

- Here is the definition of fractal dimension:
http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Shape%20Metrics/Metrics/P9%20-%20FRAC.htm

- Here is some literature:
"Fractal Analysis of the Urbanization at the Outskirts of the City: Models,
Measurement and Explanation"
http://cybergeo.revues.org/8902

- And here is the code. Don't hesitate to give advice on improvements of the
code! Maybe it would be worth developing a package for fragstats metrics?
But I still consider myself as a beginner, I don't know yet how to develop a
package

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

setwd("C:/~DOSSIER_EXPLOITATION/ADL/MR-R/~FICHES/fragstats")

# 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))

##
seqce <- seq(ncell(r))

# LISTS
area <- vector(mode="list", length=length(seqce))
perimeter <- vector(mode="list", length=length(seqce))
FRAC <- vector(mode="list", length=length(seqce))

## LOOP
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, byid=TRUE)
    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)
  }
}

perimeter <- do.call("rbind", perimeter)
area <- do.call("rbind", area)

# (P9) Fractal Dimension Index
for (i in seqce) {
  if (is.null(perimeter[[i]])) {
    FRAC[i] <- NA
  } else {
    FRAC[i] <- 2*log(0.25*perimeter[[i]])/log(area[[i]])
  }
}

FRAC <- do.call("rbind", FRAC)

## SETTING RASTER VALUES
s <- setValues(r, FRAC)

## PLOT RASTER
spplot(s)


2011/10/19 Edzer Pebesma <edzer.pebesma@uni-muenster.de>

> Sounds like something that is useful for the RStudio developers!
>
> Best wishes,
>
>
> On 10/19/2011 08:51 AM, Mathieu Rajerison wrote:
>
>> Hi Colin,
>>
>>
>> It's strange because in RStudio, there isn't any message displayed. The
>> application quits without leaving any information about the problem.
>>
>> I tried your solution and it's OK, now. I've already passed 345 windows
>> for
>> the moment.
>>
>>
>> best,
>>
>> mathieu
>>
>> 2011/10/19 Colin Rundel<rundel@gmail.com>
>>
>>  Hi Mathieu,
>>>
>>> I'm not sure why RStudio is crashing on you with that input, on my apple
>>> laptop with R 2.13.2 I get the following error when I run your sample
>>> code:
>>>
>>> Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id,
>>> "rgeos_intersection")
>>> :
>>>   TopologyException: no outgoing dirEdge found at 898287 6.27396e+06
>>>
>>> Which is how things should be handled on a topology error. It seems like
>>> there is some geos issue with treating all the buildings as a single
>>> MULTIPOLYGON. geos can choke on some geometries for any number of
>>> reasons,
>>> but my guess is that two of the buildings overlap infinitesimally after
>>> the
>>> intersection operation.
>>>
>>> My suggestion would be to use the byid option when calculating the
>>> intersection as this will find the intersection for each building
>>> separately
>>> and then wrap them back up together as a single SpatialPolygons, so for
>>> your
>>> purposes the end result should be the same.
>>>
>>> The following code worked for me:
>>>
>>> library(rgeos)
>>>
>>> load("bati.cddtes.RData")
>>> load("W.RData")
>>>
>>> z = gIntersection(bati.cddtes, W, byid=TRUE)
>>>
>>> plot(bati.cddtes)
>>> plot(W,add=TRUE)
>>> plot(z,add=TRUE,col='blue')
>>>
>>> -Colin
>>>
>>> On Oct 17, 2011, at 11:04 AM, Mathieu Rajerison wrote:
>>>
>>> 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)
>>>
>>>
>>>
>>        [[alternative HTML version deleted]]
>>
>> ______________________________**_________________
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<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<http://www.52north.org/geostatistics>
> e.pebesma@wwu.de
>

	[[alternative HTML version deleted]]

