r-lidar / lidR

Airborne LiDAR data manipulation and visualisation for forestry application
https://CRAN.R-project.org/package=lidR
GNU General Public License v3.0
587 stars 132 forks source link

"NULL value passed as symbol address" in segment_trees #622

Closed kulpojke closed 2 years ago

kulpojke commented 2 years ago

This is similar to this previous issue, which never had a full resolution, and is probably related to this, which was an issue I had previously opened, and was resolved. This is basically the same problem as my last issue, but happening with segment_trees rather than normalize_height. Like my previous issue, it works without plan(multisesson). I have attached a reproduceable example. null_val_example.zip

# import libraries
library(terra)
library(sf)
library(lidR)
library(future)

# tell lidR to use all available threads
set_lidr_threads(0)

# set up the future session
plan(multisession)

# path to data dir 
data <- "."

# path to laz files
lidar <- file.path(data, "lidar")

# path to clipping gpkg
part1 <- file.path(data, "example.gpkg")

# path to download list
downloads <- file.path(data, "downloadlist.txt")

# make dir for lidar files
dir.create(lidar, showWarnings=FALSE)

# get urls
urls <- readLines(downloads)

# download the lidar data from USGS (this seems un-necessarily complex)
for(i in seq_along(urls)){
  splt <- strsplit(urls[i], '/')
  fname <- file.path(lidar, splt[[1]][length(splt[[1]])])
  download.file(urls[i], fname, mode="wb")
}

# read the points into ctg
ctg <- readLAScatalog(lidar,
                      select="xyzcrnw",
                      filter="-drop_withheld -drop_class 9 18")

# plot
plot(ctg, chunk=TRUE)

#read gpkg, then buffer by 7m
aoi <- st_read(part1)

# verify
plot(aoi)

# make tempdir to hold results
opt_output_files(ctg) <- paste0(tempdir(), "/{XLEFT}_{YTOP}")

# clip to aoi
ctg_aoi <- clip_roi(ctg, aoi)

# rechunk, redundant?
opt_chunk_size(ctg_aoi) <- 400
opt_chunk_buffer(ctg_aoi) <- 10

# check it out
las_check(ctg_aoi)

# designate chm_path as output path
opt_output_files(ctg_aoi) <- paste0(tempdir(), "/{XLEFT}_{YTOP}")

# normalize, NOTE: Warnings about WKT CRS, probably not relevant?
ctg_norm <- normalize_height(ctg_aoi, tin())

# keep ttops in memory
opt_output_files(ctg_norm) <- ""

# locate trees
ttops <- locate_trees(ctg_norm, lmf(4), uniqueness="bitmerge")

# how many duplicates, none here, seemingly always some in larger datastes
print(paste(length(ttops$treeID) - length(unique(ttops$treeID)), "duplicates"))

# remove duplicates
#ttops <- st_cast(st_union(ttops), "POINT")

# designate chm_path as output path
opt_output_files(ctg_norm) <- paste0(tempdir(), "/{XLEFT}_{YTOP}")

# make chm
chm <- rasterize_canopy(ctg_norm, 0.5, pitfree(max_edge=c(0, 1.5)))

# designate outut path
opt_output_files(ctg_norm) <- paste0(tempdir(), "/{XLEFT}_{YTOP}")

# segmentation algorithm
algo <- dalponte2016(chm, ttops)

# segment HERE is the problem, always on first chunk to finish
ctg_segmented <- segment_trees(ctg_norm, algo)

The result:

Processing [===>--------------------------------------------------------------------------------]   5% (1/20) eta:  2mError:
NULL value passed as symbol address

The above code, and the files referenced within (downloads.txt and example.gpkg) are included in the zip

Jean-Romain commented 2 years ago

Protips: if your minimal reproducible example is more than 15 lines of code it is likely not an easily reproducible example. Here the data download failed for me. But I really appreciate you attempt at making something fully reproducible. Anyway I made my own MRE

library(lidR)

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAScatalog(LASfile)

f = tempfile(fileext = ".tif")
chm <- rasterize_canopy(las, 0.5, p2r(0.3))
terra::writeRaster(chm, f)
chm = terra::rast(f)

ttops <- locate_trees(las, lmf(4, 2))

future::plan(future::multisession)
opt_output_files(las) <- paste0(tempdir(), "/{XLEFT}_{YTOP}")
las <- segment_trees(las, dalponte2016(chm, ttops))
#> Error: NULL value passed as symbol address

This is an issue with terra, the same you previously reported and the same reported in the discussion your linked but this time with a reproducible exemple. I'll fix it with the same workaround.

Ben-Shamgochian commented 1 year ago

Hey I commented on the previous thread and am just commenting here again for visibility as this is the newer thread. What was the work around you ended up using here?

Jean-Romain commented 1 year ago

Casting to RasterLayer internally and open an issue in terra https://github.com/rspatial/terra/issues/801 so the problem is fixed by terra developper. When the new version of terra will be released I will remove the workaound