R lidR write tree points to shapefile

Question

I have created a tree point object using the locate_tree function. The tree point object is sfc_point and I want to write this to a shapefile. However, I am getting the following error.

How can I fix this?

Error:

Write error.
In addition: Warning message:
In CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :
  GDAL Error 6: Geometry type of `3D Point' not supported in shapefiles.  Type can be overridden with a layer creation option of SHPT=POINT/ARC/POLYGON/MULTIPOINT/POINTZ/ARCZ/POLYGONZ/MULTIPOINTZ/MULTIPATCH.

Code:

TTops_Man = locate_trees(LiDAR, manual(LiDAR.CHM_pitfree.TTops, radius = 2, button = "middle")) 

# to shapefile
st_write(TTops_Man ,
         "path",
         "ITD_TTops.shp",
         driver = "ESRI Shapefile")

Answer 1

The output of locate_trees() is a spatial vector with 3D coordinates for XYZ. The shapefile format does not support that. I think shapefile is an outdated format and you should better use e.g. geopackage (.gpkg) for multitple reasons.

Yet, the coordinates Z of the apices are recorded both as Z coordinates and as an attribute for "shapefile reasons". So, you can simply drop the Z coordinates with sf::st_zm()

library(lidR)
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, select = "xyz", filter = "-inside 481250 3812980 481300 3813030")
ttops <- locate_trees(las, lmf(ws = 5))
ttops
#> Simple feature collection with 29 features and 2 fields
#> Attribute-geometry relationship: 2 constant, 0 aggregate, 0 identity
#> Geometry type: POINT
#> Dimension:     XYZ
#> Bounding box:  xmin: 481260.8 ymin: 3812980 xmax: 481299.6 ymax: 3813011
#> Projected CRS: NAD83 / UTM zone 12N
#> First 10 features:
#>    treeID     Z                       geometry
#> 1       1 16.00  POINT Z (481294.7 3813011 16)
#> 2       2 26.95 POINT Z (481281.9 3813003 2...
#> 3       3 23.58 POINT Z (481278.4 3813002 2...
#> 4       4 19.75 POINT Z (481265.3 3812996 1...
#> 5       5 18.89 POINT Z (481261.8 3812999 1...
#> 6       6 16.25 POINT Z (481266.9 3813004 1...
#> 7       7 19.96 POINT Z (481273.5 3813000 1...
#> 8       8 23.70 POINT Z (481275.8 3813009 2...
#> 9       9 24.61 POINT Z (481278.2 3813011 2...
#> 10     10 23.28 POINT Z (481284.9 3812998 2...
ttops <- sf::st_zm(ttops)
ttops
#> Simple feature collection with 29 features and 2 fields
#> Attribute-geometry relationship: 2 constant, 0 aggregate, 0 identity
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 481260.8 ymin: 3812980 xmax: 481299.6 ymax: 3813011
#> Projected CRS: NAD83 / UTM zone 12N
#> First 10 features:
#>    treeID     Z                 geometry
#> 1       1 16.00 POINT (481294.7 3813011)
#> 2       2 26.95 POINT (481281.9 3813003)
#> 3       3 23.58 POINT (481278.4 3813002)
#> 4       4 19.75 POINT (481265.3 3812996)
#> 5       5 18.89 POINT (481261.8 3812999)
#> 6       6 16.25 POINT (481266.9 3813004)
#> 7       7 19.96 POINT (481273.5 3813000)
#> 8       8 23.70 POINT (481275.8 3813009)
#> 9       9 24.61 POINT (481278.2 3813011)
#> 10     10 23.28 POINT (481284.9 3812998)

^{Created on 2022-09-26 by the reprex package (v2.0.1)}