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
596 stars 131 forks source link

A missing CRS for LASCatalog, after projection(ctg), works for most things but not segmentation; problem might not be CRS #572

Closed fpegmbg closed 2 years ago

fpegmbg commented 2 years ago

Thank you for all the work you have put into developing lidR. I can often get things to work after reading your responses to error reports, but in this case, I cannot. Unfortunately, I also may not be able to post a working example because the .las files are too large--I tried to attach to this post, but unsure if they succeeded. I believe the problem is due to projections, but that could be untrue. [On second thought, I tested the code with a couple of .laz files that do have a defined CRS, and the same error occurs--slight difference, in that the tree list had duplicates, but even after distinct() it still did not work, and with the same error as below.]

I have a simple catalog of two tiles. When read in, the ctg looks like below, with NA for the CRS.

class       : LAScatalog (v1.2 format 1)
extent      : 6653000, 6655000, 1923000, 1924000 (xmin, xmax, ymin, ymax)
coord. ref. : NA 
area        : 2 kunits²
points      : 8.8 million points
density     : 4.4 points/units²
density      : 2.3 pulses/units²
num. files  : 2 

From the las metadata, I know the CRS should be EPSG:6422. I assign this with:

projection(ctg) <- 6422

Which results in:

class       : LAScatalog (v1.2 format 1)
extent      : 6653000, 6655000, 1923000, 1924000 (xmin, xmax, ymin, ymax)
coord. ref. : NAD83(2011) / California zone 4 (ftUS) 
area        : 2 kus-ft²
points      : 8.8 million points
density     : 4.4 points/us-ft²
density      : 2.3 pulses/us-ft²
num. files  : 2 

This looks correct, and works for many subsequent steps. All of the following steps work (it is already ground-classified, so I am skipping that here, although a re-classification also works), until the indicated failure. I would hope that the below code works if pointed to a local catalog replacing the paste0 expression.

THIS ALL WORKS:

ctg <- readLAScatalog(paste0(lidpath,"Tiles_0001_0002_testing"))
plot(ctg, map = TRUE)
print(ctg) #Seems to be lacking a CRS, CRS: NA
las_check(ctg, deep=TRUE)
#Projection needed for the catalog
#Is EPSG:6422, which is NAD83(2011) / California zone 4 (ftUS)
#BUT, epsg() does not work:
#Instead, projection(catalog_object) is the right approach:
#From: https://github.com/r-lidar/lidR/issues/405
projection(ctg) <- 6422
ctg
library(future)
plan(multisession)
rm(dtm)
opt_output_files(ctg) <- opt_output_files(ctg) <- paste0(tempdir(), "/{*}_dtm")
#DTM generation
dtmctg <- grid_terrain(ctg, algorithm = tin())
dtm <- rasterize_terrain(ctg, 1, tin(), overwrite=TRUE)
opt_output_files(ctg) <-  paste0(tempdir(), "/{*}_norm")
ctg_norm <- normalize_height(ctg, dtmctg)
#Area metrics on the normalized catalog:
opt_select(ctg_norm) <- "xyz"
opt_filter(ctg_norm) <- "-keep_first"
ctg_metrics_w2w <- pixel_metrics(ctg_norm, .stdmetrics_z, res = 10, pkg = "terra", overwrite=TRUE)
# plot(ctg_metrics_w2w$zmax, col = height.colors(25))
# plot(ctg_metrics_w2w$zmean, col = height.colors(25))
# plot(ctg_metrics_w2w$zmean/ctg_metrics_w2w$zmax, col = height.colors(25))
# canopyratio=ctg_metrics_w2w$zmean/ctg_metrics_w2w$zmax
# class(canopyratio)
#CHM and ttops
opt_output_files(ctg_norm) <- paste0(tempdir(), "/chm_{*}")
chm <- rasterize_canopy(ctg_norm, res = 0.5, include.lowest = FALSE, pitfree(thresholds = c(0, 10, 20), max_edge = c(0, 1.5)))
# chm <- rasterize_canopy(ctg_norm, 1, p2r(0.15))
plot(chm, col = height.colors(50))
#Treetops: Las (is in feet units, f specifies variable radius up to 25 ft)
opt_output_files(ctg_norm) <- ""
ttops <- locate_trees(ctg_norm, lmf(f), uniqueness = "bitmerge")
ttops
#This plot works, all the ttops are at locations that make sense based on the chm
plot(chm_cat, col = height.colors(50))
plot(ttops, add=TRUE)
#Now trying to segment the catalog, following the lidR book:
opt_output_files(ctg_norm) <- paste0(tempdir(), "/{*}_segmented")
#Modified the defaults for the algo because of the las units in feet:
algo <- dalponte2016(chm, ttops, th_tree=10, th_cr=0.65, max_cr=70)

THIS NEXT LINE FAILS:

ctg_segmented <- segment_trees(ctg_norm, algo)

Produces error: "Error: NULL value passed as symbol address" It also fails if I use the default settings for dalponte2016

I yesterday re-installed an entirely new R and all packages, so possibly the usual troubleshooting of being sure to update all packages may not apply. Here are the versions of several packages (R is version 4.1.3).

> packageVersion("rgdal")
[1] ‘1.5.30’
> packageVersion("sf")
[1] ‘1.0.7’
> packageVersion("lidR")
[1] ‘4.0.0’
> packageVersion("terra")
[1] ‘1.5.21’
> packageVersion("raster")
[1] ‘3.5.15’

Thank you for any insights about this error!

Jean-Romain commented 2 years ago

It looks like it has nothing to do with CRS but please share a minimal reproducible example. If you can't share a file at least share a minimal code. You report has ton of code that is not related to the problem. I can't disentangle what is related to the problem and what is not.

fpegmbg commented 2 years ago

Sorry for the poor question form. I just solved the problem, the terra library had not loaded. Segmentation works fine now that the library is here.