lasmergespatial changes polygon shape

[NatKatUP]

I have a polygon containing the ground projected crown shapes of a tree stand and want to give my point cloud points that belong to the trees their IDs.

It seems like my polygons change their shape within the funtion. It looks like the polygons would maybe get rasterized and resampled, at least their shape is much simpler and if i color the point cloud by ID in the RGL viewer, you see that not all points got the id (see screenshot).

I attach a las subset and my shape file. test.zip

library(lidR)
testlas <- readLAS("test_subset.las")
testPoly <- shapefile("EN18012_crownsPoly.shp")

plot(SpatialPoints(testlas@data[,1:2]), pch = 16, cex = 0.2)
plot(testPoly, add= TRUE, border = "red")

testlas <- lasmergespatial(testlas, testPoly, "ID") # assign ID's

plot(testlas, color = "ID") # look at the top view and compare with plot

[Jean-Romain]

To me it looks like a problem with floating point precision. You coordinates are very high number thus the issue is amplified. Moreover your polygons are tiny. The error is almost invisible for wide polygons but is very visible here. Let make a test.

library(lidR)

shp  <- shapefile("EN18012_crownsPoly.shp")
las  <- readLAS("test_subset.las")

# Shift
xmin <- min(las$X)
ymin <- min(las$Y)

# Shift everything to (0,0)
las$X <- las$X - xmin
las$Y <- las$Y - ymin
shp   <- shift(shp, x = -xmin, y = -ymin)

las <- lasmergespatial(las, shp, "ID")
plot(las, color = "ID")

It looks good. That is very annoying because the bug should occurs in rlas as well since both lidR and rlas rely on boost::geometry::covered_by.

At least you can use the following workaround waiting for a real fix.

lasmergespatialspdf = function(las, spdf, attr = NULL)
{
  xmin  <- min(las$X)
  ymin  <- min(las$Y)
  las$X <- las$X - xmin
  las$Y <- las$Y - ymin
  spdf  <- shift(spdf, x = -xmin, y = -ymin)
  las   <- lasmergespatial(las, spdf, attr)
  las$X <- las$X + xmin
  las$Y <- las$Y + ymin
  return(las)
}

[Jean-Romain]

I found the origin of the trouble. When parsed as WKT string a to small number of digit was used. I will fixed that within few minutes

spshp    <- shapefile("EN18012_crownsPoly.shp")
sfshp    <- sf::st_as_sf(spshp)
sptree1  <- spshp[1,]
sftree1  <- sfshp[1,]

# Reproduce lidR behavior
wkt <- sf::st_as_text(sftree1$geometry)
wktree1 <- rgeos::readWKT(wkt)

plot(sptree1)
plot(sftree1["ID"])
plot(wktree1)

# Use more digits
wkt <- sf::st_as_text(sftree1$geometry,  digits = 14)
wktree1 <- rgeos::readWKT(wkt)

plot(sptree1)
plot(sftree1["ID"])
plot(wktree1)

[NatKatUP]

Thanks, this works for me.

[Jean-Romain]

Fixed in v2.0.2

[NatKatUP]

Also lasclip produces errors: "No point found for within POLYGON".
Probably for the same reason.
It works with

lasclipspdf = function(las, spdf){
  xmin  <- min(las$X)
  ymin  <- min(las$Y)
  las$X <- las$X - xmin
  las$Y <- las$Y - ymin
  spdf  <- shift(spdf, x = -xmin, y = -ymin)
  las   <- lasclip(las, spdf)
  las$X <- las$X + xmin
  las$Y <- las$Y + ymin
  return(las)
}

[Jean-Romain]

I also fixed lasclip in the same commit. If you still encounter troubles it is another issue.