library(sf)
## Warning: package 'sf' was built under R version 4.0.4
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(rgdal)
## Warning: package 'rgdal' was built under R version 4.0.5
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.0.4
## 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/jmhp2/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/jmhp2/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/jmhp2/OneDrive/Documents/R/win-library/4.0/rgdal/proj
library(raster)
## Warning: package 'raster' was built under R version 4.0.4
library(tmap)
## Warning: package 'tmap' was built under R version 4.0.4
library(tmaptools)
## Warning: package 'tmaptools' was built under R version 4.0.4
library(spatialEco)
## Warning: package 'spatialEco' was built under R version 4.0.5
## 
## Attaching package: 'spatialEco'
## The following object is masked from 'package:raster':
## 
##     shift
library(RStoolbox)
## Warning: package 'RStoolbox' was built under R version 4.0.5

#Set up

setwd("C:/Users/jmhp2/Downloads/Raster/Raster")
NED <- raster("ele_a13.img")
NLCD <- raster("lc_a13.img")
slope <- raster("slope_a13.img")
stream <- raster("strm_r.img")

#Binary/Conservative Model

This model is used to show a true or false, in turn this makes it a good method to show habitats for a specific plant. Firstly we have to check that all of the imported data is useable. We made a grid populated with the imported “ele_a13.img” data, then making the grid blank by using ‘= NA’. This allowed us to fill the grid with the raster to point function that allowed us to set buffer values to the data.

blank <- NED
blank[] = NA
stream_points <- rasterToPoints(stream, spatial=TRUE)
stream1 <- distanceFromPoints(blank, stream_points) 
tm_shape(stream1)+
tm_raster(style= "cont")+
tm_layout(legend.outside = TRUE)

Next we must identify the specfics of the habitat that the model needs to show. These are:

NLCD_E <- NLCD == 42 
NLCD_M <- NLCD == 43
NLCD_EM <- NLCD_M + NLCD_E
NED_b <- NED > 1300
NED_bb <- tm_shape(NED_b)+
tm_raster(style="cat")+
tm_layout(main.title="Elevation")
slope_b <- slope > 10
slope_bb <- tm_shape(slope_b)+
tm_raster(style="cat")+
tm_layout(main.title="Slope")
NLCD_b <- NLCD_EM
ncld_b <- tm_shape(NLCD_b)+
tm_raster(style="cat")+
tm_layout(main.title="Landcover")
stream_b <- stream1 < 1000
stream_bb <- tm_shape(stream_b)+
tm_raster(style="cat")+
tm_layout(main.title="Stream")
tmap_arrange(NED_bb, slope_bb, ncld_b, stream_bb)

#Liberal/Weighted Overlay Model

This model is used to predict a range of suitablity rather than the ‘yes or no’ that coems from a binary model. The way this was done with our model was by weighted overlay. First we had to procure the grids of the data adn then apply the weight values of suitability.

NLCD_D <- NLCD == 41
NLCD_F <- NLCD_M + NLCD_E + NLCD_D

The specifics for this model are more generalized than the previous model. Here the speficiations were:

For these specification we had to rescale our raster data into values of 0 to 1.

NED_min <- cellStats(NED, stat=min)
NED_max <- cellStats(NED, stat=max)
NED_s <- ((NED - NED_min)/(NED_max - NED_min))
NED_ss <- tm_shape(NED_s)+
tm_raster(style="cont")+
tm_layout(main.title="Elevation")

slope_min <- cellStats(slope, stat=min)
slope_max <- cellStats(slope, stat=max)
slope_s <- ((slope - slope_min)/(slope_max - slope_min))
slope_ss <- tm_shape(slope_s)+
tm_raster(style="cont")+
tm_layout(main.title="Slope")

NLCD_min <- cellStats(NLCD_E, stat=min)
NLCD_max <- cellStats(NLCD_E, stat=max)
NLCD_s <- 1 - ((NLCD_E - NLCD_min)/(NLCD_max - NLCD_min))
NLCD_ss <- tm_shape(NLCD_s)+
tm_raster(style="cont")+
tm_layout(main.title="Airports")

stream_min <- cellStats(stream1, stat=min)
stream_max <- cellStats(stream1, stat=max)
stream_s <- 1- (stream1 - stream_min)/(stream_max - stream_min)
stream_ss <- tm_shape(stream_s)+
tm_raster(style="cont")+
tm_layout(main.title="Interstates")

tmap_arrange(NED_ss, slope_ss, NLCD_ss, stream_ss)