CityShadeMapper/R/rasterize_lidar_cover_class.R at 2067af5b57c92d81abd7a08ebf2c3b11999713ff · delphinekw/CityShadeMapper · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#' Get raster of cover classification from lidar points
#'
#' @param las A [lidR::LAScatalog-class()] object, or a character vector with
#' paths to LAS/LAZ objects.
#' @param res Resolution of the resulting raster.
#' @param fill.holes Logical. Try to fill holes in lidar point classification.
#' @param filename Character. Output filename. Note that if a file already exists
#' with that name, it will be overwritten.
#'
#' @return A SpatRaster with the classification of cover types:
#' 2 = ground (including low vegetation < 1m)
#' 4 = high vegetation (> 1m)
#' 6 = buildings
#' 9 = water
#' and NA values.
#' Note that points classified as bridges (class 17) will be reclassified as ground.
#' @export
#'
#' @examples
#' \dontrun{
#' pza <- system.file("extdata", "PlazaNueva.laz", package = "CityShadeMapper")
#' pza.cover <- rasterize_lidar_cover_class(pza)
#' }
rasterize_lidar_cover_class <- function(las = NULL,
                                        res = 1,
                                        fill.holes = TRUE,
                                        filename = NULL,
                                        correct.low.veg = TRUE) {

  pts <- lidR::readLAS(las, select = "xyc", filter = "-keep_class 2 3 4 5 6 9 17")
  #table(pts$Classification)
  # see classes https://github.com/Pakillo/CityShadeMapper/issues/2#issuecomment-1117648511
  # reclassify bridges as ground
  pts$Classification[pts$Classification == 17] <- as.integer(2)
  # reclassify low vegetation (<1m) as ground
  if(correct.low.veg){
    pts$Classification[pts$Classification == 3] <- as.integer(2)
  }
  # join classes 4 & 5 (vegetation > 1m high)
  pts$Classification[pts$Classification == 5] <- as.integer(4)

  pts.class <- lidR::pixel_metrics(pts,
                                   func = CityShadeMapper:::max.class(Classification),
                                   res = res)

  if (!is.null(filename)) {
    terra::writeRaster(pts.class, filename = filename, overwrite = TRUE, datatype = "INT1U")
    pts.class <- terra::rast(filename)
  }

  if (isTRUE(fill.holes)) {
    pts.class <- fill_holes(pts.class)
    if (!is.null(filename)) {
      terra::writeRaster(pts.class, filename = filename, overwrite = TRUE, datatype = "INT1U")
      pts.class <- terra::rast(filename)
    }
  }

  pts.class

}


max.class <- function(x) {max(x, na.rm = TRUE)}

fill_holes <- function(ras = NULL) {

  # try to fill holes (NA) in lidar classification of points
  ras.agg <- terra::aggregate(ras, fact = 2, fun = "modal", na.rm = TRUE)
  ras.disagg <- terra::disagg(ras.agg, fact = 2, method = "near")

}