Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
fpegmbg opened this issue Apr 19, 2022 · 2 comments

Comments

@fpegmbg
Copy link

fpegmbg commented Apr 19, 2022

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
Copy link
Collaborator

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
Copy link
Author

fpegmbg commented Apr 20, 2022

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.

@fpegmbg fpegmbg closed this as completed Apr 20, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants