R Notebook

Loading the Necessary Libraries

library(raster)

## Warning: package 'raster' was built under R version 4.0.5

## Loading required package: sp

## Warning: package 'sp' was built under R version 4.0.5

library(rgdal)

## Warning: package 'rgdal' was built under R version 4.0.5

## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/manis/OneDrive/Documents/R/win-library/4.0/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/manis/OneDrive/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/manis/OneDrive/Documents/R/win-library/4.0/rgdal/proj

library(reshape)

## Warning: package 'reshape' was built under R version 4.0.5

library(scales)

Reading the DEM (Digital Elevation Model) Data

lidar_dem <-  raster(x = "E:/Spatial analysis using R/lidar_data/BLDR_LeeHill/pre-flood/lidar/pre_DTM.tif")

#Plot the Data
plot(lidar_dem, main = "Digital Elevation Model")

Getting the Co-ordinate Reference System of the Lidar Data

crs(lidar_dem)

## CRS arguments:
##  +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs

mycrs <- crs(lidar_dem)

#X - resolution
xres(lidar_dem)

## [1] 1

#Y - resolution
yres(lidar_dem)

## [1] 1

Plot Histograms of Raster Values

Histogram represents the distribution of pixel elevation values in your data. The plot might be useful for identifying the outlier data values, assess the miin and max values in our data and explore the general distribution of elevation values of our data.

hist(lidar_dem,
     breaks  = c(1600,1800,2000,2100),
     main = "Surface Values Distribution",
     xlab = "Elevation(meters)", ylab = "Frequency",
     col = "springgreen")

## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 1%
## of the raster cells were used. 100000 values used.

Canopy height Model - It is basically the heights of the objects above the ground for this model we need DTM and DSM. We’ll load the DSM (Digital surface Model) and then try to calculate the Canopy height Model.

#load the DSM
lidar_dsm <- raster(x = "E:/Spatial analysis using R/lidar_data/BLDR_LeeHill/pre-flood/lidar/pre_DSM.tif")

#SO the formula for identifying the Canopy height Model is DSM - DTM
lidar_CHM <- lidar_dsm - lidar_dem

Plot the Canopy Height Model

#plot the Canopy height model
plot(lidar_CHM,
     breaks = c(0,2,10,20,30),
     main = "Lidar Canopy Height Model",
     col = c("white","maroon","springgreen", "darkgreen")
)

Saving the Raster File

dir.exists("data/outputs")

## [1] FALSE

if (!dir.exists("data/outputs")) {
  dir.create("data/outputs", recursive = TRUE)
}


#save the canopy height model

writeRaster(lidar_CHM,"data/outputs/lidar_chm.tiff",
            format = "GTiff",
            overwrite = TRUE)

We would need to Reclassify the Dataset in order to understand the CHM better. Raster Classification Steps are:

Data Import/cleanup
Data Exploration
More data processing and analysis
Final data analysis
Presentation

Loading and Plotting the Histogram of the Dataset

lidar_chm<- raster("E:/Spatial analysis using R/lidar_data/data/outputs/lidar_chm.tif")

hist(lidar_chm,
     main = "Distribution of raster cell values in the DTM",
     xlab = "Height (m)", #We can also use xlim wo limit our interest with respect to data
     ylab = "Number of Pixels",
     col = "red")

## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 1%
## of the raster cells were used. 100000 values used.

As we can see that as the Height Increased the Number of Pixels kept on Decreasing.

Mapping raster to new values and to reclassify we need to create a reclassification matrix. New values would be as follows:

No trees(0-2m tall) - NA
Short trees (2m - 4m) - 1
Moderate trees (4m-7m) - 2
Long trees (>7m) - 3

reclass_df <- c(0,2, NA,
                2,4, 1,
                4,7, 2,
                7, Inf, 3)


# Now we will create a matrix by specifying the columns
reclass_m <- matrix(reclass_df,
                    ncol = 3,
                    byrow = TRUE)


# Reclassify our data.
chm_classfied <- reclassify(lidar_chm,reclass_m)

barplot(chm_classfied, 
        main = "Number of pixels in each values",
        xlab = "Classes",
        ylab = "No of Pixels")

## Warning in .local(height, ...): a sample of 12.5% of the raster cells were used
## to estimate frequencies

#classify the 0 values to NA
chm_classfied[chm_classfied==0] <- NA

plot(chm_classfied,
     col = c("red","blue","green"))

#adding the custom labels
plot(chm_classfied,
     legend = FALSE,
     col = c("red","blue","green"),
     axes = FALSE,
     main = "Classifed Canopy Height Model \n Short, Medium, Tall trees")

chm_colors <- c("red","blue","green")

legend("topright",
  legend = c("short trees", "moderate trees", "tall trees"),
  fill = chm_colors,
  border = FALSE,
  bty = "n" #turns off the legend border
  )