Tarea 3 LiDAR-ALS
Luis Castillo
2022-09-05
Cree un script (código) que resuelva y muestre las siguientes preguntas (cada pregunta otorga 5 puntos a su tarea de calificación):
- Lea el archivo 001134.las en R. Subconjunto del data.frame de este objeto LAS y muestre la estructura, str(), de estos datos. Proyecte este LAS en UTM para la zona 18 al norte y luego genere un LAS normalizado. Trazar la altura y la intensidad utilizando el gráfico del paquete “lidR”. Genere el modelo de terreno digital y el modelo de elevación digital de la parte superior del dosel. Exporte este ráster fuera de R; necesita usar la función writeRaster() del ráster del paquete. Debe revisar el uso de la función writeRaster() utilizando la función help(); ayuda (escribir trama).
# To clean environment
rm(list = ls(all.names = TRUE))
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 518315 27.7 1158800 61.9 644200 34.5
## Vcells 924256 7.1 8388608 64.0 1634527 12.5list.of.packages <- c("rgdal", "rgeos", "raster", "sp", "lidR", "rgl")# To create a list of the packages.
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)>0) install.packages(new.packages) #To install the absent packageslibrary(rgdal)## Loading required package: sp## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/luisc/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/luisc/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library(rgeos)## rgeos version: 0.5-9, (SVN revision 684)
## GEOS runtime version: 3.9.1-CAPI-1.14.2
## Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.5-0
## Polygon checking: TRUElibrary(raster)
library(sp)
library(lidR )##
## Attaching package: 'lidR'## The following objects are masked from 'package:raster':
##
## projection, projection<-Read the file 001134.las in R. Subset the data.frame of this LAS object and show the structure, str(), of this data.
lidar_ALS_data_1 = lidR::readLAS(files = "C:/Users/luisc/Downloads/001134.las")
print(lidar_ALS_data_1)#To see the general information of the.laz## class : LAS (v1.2 format 1)
## memory : 174.1 Mb
## extent : 277000, 278000, 567000, 567480.5 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## area : 484000 units²
## points : 3.8 million points
## density : 7.86 points/units²
## density : 5.49 pulses/units²str(lidar_ALS_data_1@data)## Classes 'data.table' and 'data.frame': 3802802 obs. of 16 variables:
## $ X : num 278000 278000 278000 278000 278000 ...
## $ Y : num 567057 567060 567059 567055 567058 ...
## $ Z : num 24 31.9 31.3 23.6 30.9 ...
## $ gpstime : num 247890 247890 247890 247890 247890 ...
## $ Intensity : int 1 0 1 0 1 0 2 2 0 2 ...
## $ ReturnNumber : int 1 1 1 2 1 2 1 1 2 1 ...
## $ NumberOfReturns : int 1 1 2 2 2 2 1 1 2 1 ...
## $ ScanDirectionFlag: int 0 0 0 0 0 0 0 0 0 0 ...
## $ EdgeOfFlightline : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Classification : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Synthetic_flag : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ Keypoint_flag : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ Withheld_flag : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ ScanAngleRank : int 0 0 0 0 0 0 0 0 0 0 ...
## $ UserData : int 0 0 0 0 0 0 0 0 0 0 ...
## $ PointSourceID : int 0 0 0 0 0 0 0 0 0 0 ...
## - attr(*, ".internal.selfref")=<externalptr>Project this LAS in UTM for the zone 18 at north
projStr = "+proj=utm +zone=18 +north +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs + datum=WGS84"
crs(lidar_ALS_data_1) <- projStr## Warning: EPSG code not found: header not updated. Try to provide the ESPG code
## instead of a characterlidar_ALS_data_1@crs## Coordinate Reference System:
## User input: +proj=utm +zone=18 +north +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs + datum=WGS84
## wkt:
## BOUNDCRS[
## SOURCECRS[
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["UTM zone 18N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-75,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",16018]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Position Vector transformation (geog2D domain)",
## ID["EPSG",9606]],
## PARAMETER["X-axis translation",0,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",0,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",0,
## ID["EPSG",8607]],
## PARAMETER["X-axis rotation",0,
## ID["EPSG",8608]],
## PARAMETER["Y-axis rotation",0,
## ID["EPSG",8609]],
## PARAMETER["Z-axis rotation",0,
## ID["EPSG",8610]],
## PARAMETER["Scale difference",1,
## ID["EPSG",8611]]]]the generate a normalize LAS
lidar_ALS_dtm = lidR::grid_terrain(las=lidar_ALS_data_1, res = 10, algorithm = knnidw(k = 6, p = 2))plot(lidar_ALS_dtm, main = 'Modelo digital de terreno (DTM)')
lidar_ALS_data_1_normal = normalize_height(las=lidar_ALS_data_1, lidar_ALS_dtm)lidR::plot(lidar_ALS_data_1_normal)Generate the digital terrain model and the digital elevation model of the top of the canopy
lidar_ALS_CH = lidR::grid_canopy(las=lidar_ALS_data_1_normal, res = 10, p2r(subcircle = 0.2, na.fill = knnidw(k = 6, p = 2)))
plot(lidar_ALS_CH, main = 'Canopy Height Model (CHM)')
Plot the height and Intensity using the plot of the package “lidR”.
lidR::plot(lidar_ALS_data_1_normal, color = "Intensity", breaks = 'quantile', nbreaks = 5) #Intesidad
lidR::plot(lidar_ALS_data_1_normal) #AlturaExport this raster out of R; you need to use the function writeRaster() of the package raster. You need to review the use of the function writeRaster() using the function help(); help(writeRaster).
help("writeRaster")## starting httpd help server ... donewriteRaster(lidar_ALS_CH,"lidar.tif", overwrite=TRUE)- Lea el archivo 001132.las en R. Proyecte este LAS en UTM para la zona 18 al norte y luego genere un LAS normalizado. Elimina el ruido usando la función filter_noise_LAS(), creamos esta función en el tutorial. Genere el ráster de las siguientes dos métricas usando la función grid_metrics(): 1) el ráster del máximo de la altura y 2) el ráster de la mediana de la altura.
Read the file 001132.las in R. Project this LAS in UTM for the zone 18 at north
LIDAR_2 = lidR::readLAS(files = "C:/Users/luisc/Downloads/001132.las")projStr = "+proj=utm +zone=18 +north +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs + datum=WGS84"
crs(LIDAR_2) <- projStr## Warning: EPSG code not found: header not updated. Try to provide the ESPG code
## instead of a characterLIDAR_2@crs## Coordinate Reference System:
## User input: +proj=utm +zone=18 +north +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs + datum=WGS84
## wkt:
## BOUNDCRS[
## SOURCECRS[
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["UTM zone 18N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-75,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",16018]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Position Vector transformation (geog2D domain)",
## ID["EPSG",9606]],
## PARAMETER["X-axis translation",0,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",0,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",0,
## ID["EPSG",8607]],
## PARAMETER["X-axis rotation",0,
## ID["EPSG",8608]],
## PARAMETER["Y-axis rotation",0,
## ID["EPSG",8609]],
## PARAMETER["Z-axis rotation",0,
## ID["EPSG",8610]],
## PARAMETER["Scale difference",1,
## ID["EPSG",8611]]]]normalize LAS
lidar2_DTM = lidR::grid_terrain(las=LIDAR_2, res = 10, algorithm = knnidw(k = 6, p = 2))
plot(lidar2_DTM, main = 'Modelo Digital de Terreno (DTM)')
lidar_ALS_Norm2 = normalize_height(LIDAR_2, lidar2_DTM)Eliminate the noise using the function filter_noise_LAS(), we create this function in the tutorial.
filter_noise_LAS = function(las)
{
p95 <- grid_metrics(lidar_ALS_Norm2, func = ~quantile(Z, probs = 0.95), res = 10)
las = merge_spatial(las, p95, "p95")
las = filter_poi(las, Z < p95)
las$p95 = NULL
return(las)
}
IM_filtrada = filter_noise_LAS(lidar_ALS_Norm2)
lidR::plot(IM_filtrada)Generate the raster of the next two metrics using the function grid_metrics(): 1) the raster of maximum of the height and 2) the raster of median of the height.
fun_max_mean = function(x) {list(max = max(x), median = median(x))}
Z_max_mean = grid_metrics(IM_filtrada, ~ fun_max_mean (Z), res = 10)
plot(Z_max_mean )