Search code examples
rspatstatr-maptoolsr-ppp

Calculate Area of Discrete Regions in a Raster Image Using R

I have a raster image (20231031_110019.png, although it could be any format, jpg, png, bmp) and would like to individually enter image description herecalculate the areas of discrete regions (I will call them 'islands' although they are actually areas of fungal growth on a petri dish) in this raster.

I am trying to use the spatstat package in R to do this. I found a post stating the following code can do this:

library(spatstat)
## create example data
dd <- dilation(redwood, 0.5, polygonal=FALSE)
## find connected components
P <- connected(dd)
## convert to a tessellation
B <- tess(image=P)
## compute area of each tile
answer <- tile.areas(B)

An owin object (redwood in the above instance) is needed to execute the dilation command and I cannot find a way to convert a raster into an owin object without the maptools package, which is no longer available on CRAN. I also have been unable to install an archived version maptools using devtools.

What would a workaround be for this?

Is there another way to convert a raster into a useable object?

Is there an equally straightforward method for calculating the areas of the fungal 'islands' in my raster image?

Many thanks!

Emily

Solution

  • With "terra" you can first identify the patches, and then count the number of cells in each patch. You could then multiply that number with the actual size of each pixel to get areas. "patches" are defined as islands of values in a sea of missing values; so in this case 255 needs to be set as the missing values flag.

    library(terra)
x <- rast("https://i.sstatic.net/XjIVG.png")
NAflag(x) <- 255
p <- patches(x)
f <- freq(p)
f
#  layer value count
#1     1     1 36015
#2     1     2 17132
#3     1     3 15571
#4     1     4 15345
#5     1     5 38977
#6     1     6 48821
#7     1     7 49660
#8     1     8 41246

    Some illustration

    bnds <- as.polygons(p)
area <- subst(p, f$value, f$count)

plot(trim(p, 100), main="ID")
plot(trim(area, 100), main="area")
lines(bnds, col="dark gray")
text(bnds, cex=.75)

    enter image description here

    To go from pixel counts to area you need to know the size of your image in the real world. Let's say it has a height of 20 and a width of 15 cm. In that case you could do

    ext(x) <- c(0, 15, 0, 20)
cell_size <- prod(res(x))

f$area <- f$count * cell_size
head(f)
#  layer value count     area
#1     1     1 36015 0.900375
#2     1     2 17132 0.428300
#3     1     3 15571 0.389275
#4     1     4 15345 0.383625
#5     1     5 38977 0.974425
#6     1     6 48821 1.220525

    With area expressed in cm2.

    At that point you could also go the polygons-route suggested by Eduardo like this:

    library(terra)
x <- rast("https://i.sstatic.net/XjIVG.png")
NAflag(x) <- 255
ext(x) <- c(0, 15, 0, 20)
crs(x) <- "local"

a <- as.polygons(x)  |> disagg()
a$id <- 1:nrow(a)
a$area <- expanse(a, transform=FALSE)


plot(a, "area")
text(a)

head(data.frame(a))
#  XjIVG id     area
#1     0  1 0.900375
#2     0  2 0.428300
#3     0  3 0.389275
#4     0  4 0.383625
#5     0  5 0.974425
#6     0  6 1.220525

    PS: Note that if this were geographic data you would want to set the coordinate reference system with crs<- and use cellSize to get the actual size of the cells in the real world and then use zonal to compute the area of each patch.