identicalCRS(x, y) is not TRUE error in R for individual tree detection using chm

Question

I have a LiDAR point cloud data for a dimension of 250*250 m^2 area (a forest region). I need to separate out individual trees using that data.

I created Canopy Height Model(CHM) using LASTools and used that CHM for tree delineation. I'm attaching that chm file (this raster will give the height information)

I tried to use rLiDAR package available in R.

I coded like this

library(rLiDAR)
schm <- CHMsmoothing(chm, "mean", 5)

# Setting the fws:
fws <- 5 # dimention 5x5

# Setting the specified height above ground for detectionbreak
minht <- 8.0

# Getting the individual tree detection list 
treeList <- FindTreesCHM(schm, fws, minht)

But it's giving an error

Error: identicalCRS(x, y) is not TRUE

How can I overcome this?

Answer 1

In function FindTreesCHM, at lines 17-18, we find:

XYmax <- SpatialPoints(xyFromCell(setNull, Which(setNull == 
    1, cells = TRUE)))

Which creates a SpatialPoints. The problem is that object has no projection set:

projection(XYmax)
# [1] NA

Then, the line 19

htExtract <- over(XYmax, as(chm, "SpatialGridDataFrame"))

Throws an error because XYmax has no projection set, while chm has:

projection(chm)
# [1] "+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

And as function over first check the projections of objects, we get the error:

identicalCRS(XYmax, as(chm, "SpatialGridDataFrame"))
# [1] FALSE

A workaround would be to write your own function, adding a line setting the projection of XYmax to the projection of chm.

Also, there is an error thrown by the line 22 because of line 21.

This function can be easily fixed, but I would highly recommend to contact the maintainer of the package (maintainer("rLiDAR")).

---

Here is one possible fix:

library(rLiDAR)
library(raster)

FindTreesCHM.fix <- function(chm, fws = 5, minht = 1.37) 
{
    if (class(chm)[1] != "RasterLayer") {
        chm <- raster(chm)
    }
    if (class(fws) != "numeric") {
        stop("The fws parameter is invalid. It is not a numeric input")
    }
    if (class(minht) != "numeric") {
        stop("The minht parameter is invalid. It is not a numeric input")
    }
    w <- matrix(c(rep(1, fws * fws)), nrow = fws, ncol = fws)
    chm[chm < minht] <- NA
    f <- function(chm) max(chm)
    rlocalmax <- focal(chm, fun = f, w = w, pad = TRUE, padValue = NA)
    setNull <- chm == rlocalmax
    XYmax <- SpatialPoints(xyFromCell(setNull, Which(setNull == 
        1, cells = TRUE)))
    projection(XYmax) <- projection(chm)
    htExtract <- over(XYmax, as(chm, "SpatialGridDataFrame"))
    treeList <- cbind(slot(XYmax, "coords"), htExtract)
    colnames(treeList) <- c("x", "y", "height")
    return(treeList)
}


chm <- raster("dem_test.tif")
schm <- CHMsmoothing(chm, "mean", 5)
fws <- 5
minht <- 8.0
treeList <- FindTreesCHM.fix(schm, fws, minht)

#           x       y  height
# 1  256886.5 4110940 14.1200
# 2  256805.5 4110884 13.8384
# 3  256756.5 4110880 15.2004
# 4  256735.5 4110874 17.6100
# 5  256747.5 4110840 18.2592
# 6  256755.5 4110828 19.9252
# 7  256777.5 4110806 12.7180
# 8  256780.5 4110802 14.6512
# 9  256780.5 4110792 15.8532
# 10 256763.5 4110786 18.7128
# 11 256766.5 4110764 14.4972