CODEWITHSUNDEEP
rspatstatr-maptoolsr-ppp
Emily Rowlands
Emily Rowlands

Reputation: 31

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

Upvotes: 0

Views: 452

Answers (6)

Robert Hijmans
Robert Hijmans

Reputation: 47016

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.

Upvotes: 4

Ege Rubak
Ege Rubak

Reputation: 4507

You don’t need maptools since your data is in png format. The following modification of the spatstat workflow should work (png package must be installed):

library(spatstat)
# Read in png as a matrix:
mat_png <- png::readPNG("XjIVG.png")
# There are different conventions of how a matrix should be layed out to 
# represent points in a Cartesian coordinate system, so a conversion is needed:
mat_im <- transmat(mat_png, from = "European", to = "spatstat")
i <- as.im(mat_im)
# White pixels are 1 and black 0, so 1 should be treated as background:
cc <- connected.im(i, background = 1)
B <- tess(image=cc)
answer <- tile.areas(B)
# Areas of the components in number of pixels (multiply by pixel area to get
# physical area):
answer
#>     1     2     3     4     5     6     7     8 
#> 41246 49660 48821 38977 15345 15571 36015 17132
# Visual confirmation of components:
plot(B)

Upvotes: 3

datawookie
datawookie

Reputation: 6524

If you want to use {spatstat} then try this. I couldn't get it to work with owin but a point pattern seemed feasible.

IMG <- "XjIVG-smaller.png"

library(spatstat)
library(png)

image <- readPNG(IMG)

coords <- which(image == 1, arr.ind = TRUE)
coords <- as.data.frame(coords)
points <- ppp(coords$row, coords$col, window=owin(range(coords$row), range(coords$col)))

radius <- 5
dilated_points <- dilation(points, radius)

plot(dilated_points)

Upvotes: 0

datawookie
datawookie

Reputation: 6524

Try using the {mmand} package.

IMG <- "XjIVG.png"

library(png)
library(mmand)

image <- readPNG(IMG)

display(image)

kernel <- shapeKernel(c(10, 10), type="box")

erode(image, kernel)

dilate(image, kernel)

The eroded image:

The eroded image.

The dilated image:

The dilated image.

🚨 The operations are being applied to the white areas of the image. If you want to erode/dilate the black bits then you need to either negate the image or just use the "other" operation (dilation on black is like erosion on white).

Upvotes: 0

Eduardo Moreira
Eduardo Moreira

Reputation: 101

I suggest a different approach. Instead of 'spatstat', consider using the 'terra' package to read the raster and vectorise it, with the functions 'rast', and as 'as.polygons'. Then you can convert to an 'sf' object, with the function 'st_as_sf', and use the function 'st_cast', which will convert the multipart polygon into multiple singlepart polygons. Finally, you can use the function 'st_area' to calculate the area of all "islands". Please see the example code below.

library(terra)
library(sf)

# Read the raster into a SpatRaster object
raster_i <- rast("XjIVG.png")

# Chek the raster resolution
# It by defaut shoud be 1x1
res(raster_i)

# Polygonyze the raster
poly_i <- as.polygons(raster_i, aggregate=T)

# Separate the polygons into individual features
poly_i <- st_cast(st_as_sf(poly_i[1]), "POLYGON")

# Calculate the features area
# Since the resolution is 1x1 the area will be in pixel squared. 
# You can multiply by your scale factor squared to convert the pixel squared units to metric units.
poly_i$area <- st_area(poly_i)

# write the area values of all "islands" in a .csv table
write.csv(data.frame(poly_i)[,c(1,3)], "poly_i.csv")

# Plot the features with a color scale by area
plot(poly_i["area"])

The result of this code will be a 'sf' vector object with the islands' area in pixel-squared units. To convert to your measurement units, you can multiply these values by your scale factor squared. See the results plotted with a colour scale by the "islands" area: enter image description here

If this was helpful, please remember to recommend the answer.

Upvotes: 5

IRTFM
IRTFM

Reputation: 263301

You can still download and install the {maptools} package from the CRAN archive.

> install.packages("https://cran.r-project.org/src/contrib/Archive/maptools/maptools_1.1-8.tar.gz", repo=NULL, type="source")
> library(maptools)

Upvotes: 1

Related Questions