Problem with `lsm_l_ai()`

Today, I unexpectedly came across a problem with lsm_l_ai() - it gave NA on some of my small raster.
My raster file - test.zip.

library(raster)
library(landscapemetrics)

x = raster("test.tif")
lsm_l_ai(x)

#> # A tibble: 1 x 6
#>   layer level     class    id metric value
#>   <int> <chr>     <int> <int> <chr>  <dbl>
#> 1     1 landscape    NA    NA ai        NA

landscapemetrics/R/lsm_c_ai.R

Line 163 in d87ebd2

ai[is.nan(ai)] <- NA

Shouldn't the aggregation be 100 for only one-cell classes? @mhesselbarth @marcosci

Edit:
"The maximum aggregation is achieved when the patch type consists of a single, compact patch, which is not necessarily a square patch."

Edit2:
"AI is undefined and reported as “N/A" in the "basename".class file if the class consists of a single cell. "

Ok, based on the last quote, it seams that the landscapemetrics implementation is correct. However, I am still unsure if this approach is the best.

Edit:
I've learned that this metric is very dependable on the input data resolution:

library(raster)
#> Loading required package: sp
library(landscapemetrics)

x = raster("test2.tif")
x2 = disaggregate(x, 2)
x2
#> class       : RasterLayer 
#> dimensions  : 60, 60, 3600  (nrow, ncol, ncell)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : 0, 30, 0, 30  (xmin, xmax, ymin, ymax)
#> coord. ref. : NA 
#> data source : in memory
#> names       : test2 
#> values      : 1, 2  (min, max)
x3 = disaggregate(x, 8)
x3
#> class       : RasterLayer 
#> dimensions  : 240, 240, 57600  (nrow, ncol, ncell)
#> resolution  : 0.125, 0.125  (x, y)
#> extent      : 0, 30, 0, 30  (xmin, xmax, ymin, ymax)
#> coord. ref. : NA 
#> data source : in memory
#> names       : test2 
#> values      : 1, 2  (min, max)

lsm_c_ai(x)
#> # A tibble: 2 x 6
#>   layer level class    id metric value
#>   <int> <chr> <int> <int> <chr>  <dbl>
#> 1     1 class     1    NA ai      97.1
#> 2     1 class     2    NA ai      94.3
lsm_c_ai(x2)
#> # A tibble: 2 x 6
#>   layer level class    id metric value
#>   <int> <chr> <int> <int> <chr>  <dbl>
#> 1     1 class     1    NA ai      98.6
#> 2     1 class     2    NA ai      97.2
lsm_c_ai(x3)
#> # A tibble: 2 x 6
#>   layer level class    id metric value
#>   <int> <chr> <int> <int> <chr>  <dbl>
#> 1     1 class     1    NA ai      99.6
#> 2     1 class     2    NA ai      99.3


y = raster("test.tif")
y
#> class       : RasterLayer 
#> dimensions  : 30, 30, 900  (nrow, ncol, ncell)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 30, 0, 30  (xmin, xmax, ymin, ymax)
#> coord. ref. : NA 
#> data source : in memory
#> names       : test 
#> values      : 1, 4  (min, max)
y2 = disaggregate(y, 2)
y2
#> class       : RasterLayer 
#> dimensions  : 60, 60, 3600  (nrow, ncol, ncell)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : 0, 30, 0, 30  (xmin, xmax, ymin, ymax)
#> coord. ref. : NA 
#> data source : in memory
#> names       : test 
#> values      : 1, 4  (min, max)

lsm_c_ai(y)
#> # A tibble: 3 x 6
#>   layer level class    id metric value
#>   <int> <chr> <int> <int> <chr>  <dbl>
#> 1     1 class     1    NA ai      92.7
#> 2     1 class     2    NA ai      94.3
#> 3     1 class     4    NA ai      NA
lsm_c_ai(y2)
#> # A tibble: 3 x 6
#>   layer level class    id metric value
#>   <int> <chr> <int> <int> <chr>  <dbl>
#> 1     1 class     1    NA ai      96.4
#> 2     1 class     2    NA ai      97.2
#> 3     1 class     4    NA ai     100

^{Created on 2019-02-24 by the reprex package (v0.2.1)}

test2.zip
test.zip

Hi @mhesselbarth @bitbacchus,
this issue bit me again. I am calculating landscape metrics for hundreds of landscapes. Some of them have one-cell categories, and thus I get NA for lsm_l_ai(). Could one of you try to calculate AI on the landscape level based on the data in the following example using fragstat and let me know if the result is also NA?

library(raster)
#> Loading required package: sp
library(landscapemetrics)
data("landscape")

# original data -----------------------------------------------------------
plot(landscape, col = c("red", "green", "blue"))

lsm_l_ai(landscape)
#> # A tibble: 1 x 6
#>   layer level     class    id metric value
#>   <int> <chr>     <int> <int> <chr>  <dbl>
#> 1     1 landscape    NA    NA ai      81.1

# data with a one-cell class ----------------------------------------------
landscape[444] = 4

plot(landscape, col = c("red", "green", "blue", "yellow"))

lsm_l_ai(landscape)
#> # A tibble: 1 x 6
#>   layer level     class    id metric value
#>   <int> <chr>     <int> <int> <chr>  <dbl>
#> 1     1 landscape    NA    NA ai        NA

^{Created on 2020-12-31 by the reprex package (v0.3.0)}

I don't have a Windows machine available anymore, sorry.

Reading the documentation again, it kind of makes sense that it is NA for only one cell though since there are no adjacencies at all in this case, which is different to a fully aggregated patch of several cells.

Oh...but I see the problem now. It actually is a strange behavior if more than 1 class exists as in your example above. I think the problem is that we need to the take the sum of AI on class level, which will be NA if one class (in your example class4) is NA. We could simply add na.rm = TRUE to the sum. Would be indeed interesting what FRAGSTATS returns...

library(landscapemetrics)

landscape[444] = 4

lsm_l_ai(landscape)
#> # A tibble: 1 x 6
#>   layer level     class    id metric value
#>   <int> <chr>     <int> <int> <chr>  <dbl>
#> 1     1 landscape    NA    NA ai        NA

# get aggregation index for each class
ai <- lsm_c_ai(landscape)

# get proportional class area
pland <- lsm_c_pland(landscape,
                     directions = 8)

# final AI index
(sum(ai$value * (pland$value / 100), na.rm = TRUE))
#> [1] 80.82197

^{Created on 2021-01-05 by the reprex package (v0.3.0)}

@mhesselbarth I just remembered that I have a remote access to some windows computers at my Uni.

I created a tiny simple raster:

library(landscapemetrics)
library(raster)
#> Loading required package: sp
r = raster(matrix(c(rep(1, 8), 2, rep(1, 7)), ncol = 4))

lsm_l_ai(r)
#> # A tibble: 1 x 6
#>   layer level     class    id metric value
#>   <int> <chr>     <int> <int> <chr>  <dbl>
#> 1     1 landscape    NA    NA ai        NA

It returns NA in landscapemetrics, however, when I calculated AI in Fragstat, the result was:

This is consistent with your idea above:

library(landscapemetrics)
library(raster)
#> Loading required package: sp
r <- raster(matrix(c(rep(1, 8), 2, rep(1, 7)), ncol = 4))

# get aggregation index for each class
ai <- lsm_c_ai(r)

# get proportional class area
pland <- lsm_c_pland(r, directions = 8)

# final AI index
(sum(ai$value * (pland$value / 100), na.rm = TRUE))
#> [1] 89.48864

^{Created on 2021-01-06 by the reprex package (v0.3.0)}

Okay, so that's an easy fix. I can do that later today.

Can you please create a raster with three classes (2 classes with more than one cell and 1 class with only one cell) and calculate AI on class level as well just to make sure the results will be numeric for two classes and NA for one class

library(raster)
#> Loading required package: sp
r2 = raster(matrix(c(1, 1, 1, 3, 3, 3, 3, 3, 2, 1, 1, 3, 3, 1, 1, 1), ncol = 4))
plot(r2)

Fragstats results:

^{Created on 2021-01-06 by the reprex package (v0.3.0)}

017202b
Perfect. Should be fixed!

library(landscapemetrics)
library(raster)

r2 = raster(matrix(c(1, 1, 1, 3, 3, 3, 3, 
                     3, 2, 1, 1, 3, 3, 1, 1, 1), ncol = 4))

lsm_c_ai(r2)
#> # A tibble: 3 x 6
#>   layer level class    id metric value
#>   <int> <chr> <int> <int> <chr>  <dbl>
#> 1     1 class     1    NA ai      70  
#> 2     1 class     2    NA ai      NA  
#> 3     1 class     3    NA ai      62.5
lsm_l_ai(r2)
#> # A tibble: 1 x 6
#>   layer level     class    id metric value
#>   <int> <chr>     <int> <int> <chr>  <dbl>
#> 1     1 landscape    NA    NA ai      62.3

^{Created on 2021-01-06 by the reprex package (v0.3.0)}

Awesome. Thanks Max!