LinearDistScript.Rmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---

## Preparing Linear Measures

First we will begin with determining the distance between several features. In our first example, we want to measure distance from each mule deer location to the nearest stream if it is determined a priori that water or riparian habitats influence mule deer distribution in our study area. While this may not seem like a very complicated process, there are numerous steps needed to achieve this feat. We will need to use the package spatstat that will help us in creating individual segments with nodes for linear features such as roads and streams/rivers.

1\. Open the script LinearDistscript.Rmd" and run code directly from the script

2\. First we need to load the packages needed for the exercise

```{r warning=FALSE, message=FALSE}
library(spatstat)
library(polyCub)#replaces gpclib function
library(OneR)#bin function
library(sf)
library(terra)
library(dplyr)
library(nngeo)
library(rlang)
```

3\. Now let's have a separate section of code to include projection information we will use throughout the exercise. In previous versions, these lines of code were within each block of code

```{r warning=FALSE, message=FALSE}
ll.crs <- st_crs(4269)
utm.crs <- st_crs(9001)
albers.crs <- st_crs(5070)
```

4\. Load the mule deer dataset we used in the previous exercise

```{r}
muleys <-read.csv("data/muleysexample.csv", header=T)
#Remove outlier locations
coords <- st_as_sf(muleys, coords = c("Long", "Lat"), crs = ll.crs)
plot(st_geometry(coords),axes=T)
deer.spdf <- st_crop(coords, xmin=-107.0,xmax=-110.5,ymin=37.8,ymax=39.0)#Visually identified based on previous plot
plot(st_geometry(deer.spdf),axes=T)

#Project deer.spdf to Albers as in previous exercise
deer.albers <-st_transform(deer.spdf, crs=albers.crs)
plot(st_geometry(deer.albers,axes=T))


```

5.  Only use code in next section for example exercise so fewer locations are used.

```{r}
muleys <- muleys[sample(nrow(muleys), 100),]

#Make a spatial data frame of locations after removing outliers
coords <- st_as_sf(muleys, coords = c("Long", "Lat"), crs = ll.crs)
plot(st_geometry(coords),axes=T)
deer.spdf <- st_crop(coords, xmin=-107.0,xmax=-110.5,ymin=37.8,ymax=39.0)#Visually identified based on previous 
#Project deer.spdf to Albers as in previous exercise
deer.albers <-st_transform(deer.spdf, crs=albers.crs)
plot(st_geometry(deer.albers),axes=T)
```

6\. If we get some NA errors because our grid does not encompass our panther locations then we can expand the grid size extending beyond our locations using methods in an earlier exercise.

```{r}
# Create vectors of the x and y points using boundary box created around deer locations
bb1 <- st_bbox(deer.albers)
     
increment = 1000
minx=(min(bb1$xmin)-(increment))
maxx=(max(bb1$xmax)+(increment))
miny=(min(bb1$ymin)-(increment))
maxy=(max(bb1$ymax)+(increment))

my_bbox = st_bbox(c(xmin = minx, xmax = maxx, 
                    ymin = miny, ymax = maxy),
                  crs = 5070)

AlbersSP <- st_as_sfc(my_bbox)
```

7\. Load the necessary road and rivers shapefiles already in Albers projection to match previous vegetation raster.

```{r message=FALSE,warning=FALSE}
roads<-st_read("data/AlbersRoads.shp")
rivers<-st_read("data/AlbersRivers.shp")
```

10\. Load vegetation raster layer tif that came in the Albers projection from the online source.

```{r}
veg <-rast("data/cropnlcd.tif")
```

```{r}
#Check to see all our layers are now in Albers projection
st_crs(veg)
st_crs(deer.albers)
st_crs(AlbersSP)
```

```{r}
plot(veg)
plot(st_geometry(deer.albers),add=T, col="red")
```

11\. Then we need to expand the bounding polygon so all locations are included. We can then make the bounding polygon (AlbersSP) a class owin in order to proceed with functions in package spatstat.

```{r warning=FALSE,message=FALSE}
buffSP <- st_buffer(AlbersSP,2000)
plot(st_geometry(buffSP))
plot(st_geometry(deer.albers),add=T,col="red")
```

12\. Code below will be for use with the spatstat package to convert segments of line layers (e.g., roads, rivers) to lines to enable distance to feature from deer locations. Most calculations with spatstat require 3 new classes so most code is created to achieve this goal:

"owin" Observation windows "ppp" Planar point patterns "psp" Planar segment patterns

```{r}
#Replace AlbersSP with buffSP if using a subsample of deer
bdy.owin <- as.owin(buffSP)
is.owin(bdy.owin)
#It is TRUE so now we can move forward with the analysis
```

12\. Now clip the raster using the buffered bounding box (buffSP) created in step 5.

```{r warning=FALSE,message=FALSE}
cropSP <- st_buffer(AlbersSP,1000)
bbclip <- crop(veg, cropSP) 
cliproads <- st_intersection(roads, cropSP, byid=TRUE)
cliprivers <- st_intersection(rivers, cropSP, byid=TRUE)

plot(st_geometry(buffSP))
plot(bbclip,add=T)
plot(st_geometry(cropSP),add=T)
plot(st_geometry(deer.albers),add=T,col="red")
plot(st_geometry(cliproads),add=T)
plot(st_geometry(cliprivers), col="blue",add=T)
```

13\. We will start with the road layer by converting a single line to a set of segments packaged as a function.

```{r warning=FALSE,message=FALSE, eval=FALSE}
segs <- st_segments(cliproads,units::set_units(50,m))
segs_lines <- segs %>%
  st_as_sf(coords = c("x", "y"), crs(deer.albers)) %>%
  group_by(LINEARID) %>%
  dplyr::summarize() %>%  
  filter(st_is_empty(.) == FALSE) %>%
  st_cast("MULTILINESTRING") 

plot(st_geometry(segs_lines))
segs.psp <- as.psp(segs_lines, window=bdy.owin)
#The segments as a planar segment pattern:
plot(segs.psp)
points(deer.albers)
segs.psp[1:5]
```

14\. We first need to handle the mule deer locations. We need to make mule deer xy coordinates a planar point pattern (i.e., ppp) for use in package spatstat.

```{r warning=FALSE,message=FALSE, eval=FALSE}
newdeer <- data.frame(st_coordinates(deer.albers))
xy.ppp <- as.ppp(newdeer,W=bdy.owin)
plot(xy.ppp)
```

15\. Now we can determine the distance from mule deer locations (xy.ppp) to the nearest road

```{r eval=FALSE}
roaddist <- nncross(xy.ppp, segs.psp)$dist
#Or identify segment number closest to each point
v <- nearestsegment(xy.ppp,segs.psp)#Identifies segment number not a distance
plot(segs.psp)
plot(xy.ppp[1], add=TRUE, col="red")
plot(segs.psp[v[1]], add=TRUE, lwd=8,col="blue")
```

16\. Now we do the same to a river layer by converting a single line to a set of segments packaged as a function.

```{r eval=FALSE}
segs2 <- st_segments(cliprivers,units::set_units(50,m))
segs2_lines <- segs2 %>%
  st_as_sf(coords = c("x", "y"), crs(deer.albers)) %>%
  group_by(LINEARID) %>%
  dplyr::summarize() %>%  
  filter(st_is_empty(.) == FALSE) %>%
  st_cast("MULTILINESTRING") 

plot(st_geometry(segs2_lines))
segs2.psp <- as.psp(segs2_lines, window=bdy.owin)
#The segments as a planar segment pattern:
plot(segs2.psp)
points(deer.albers)
is.psp(segs2.psp)
#All is TRUE so now we can move forward with the analysis
```

17\. Now we can determine the distance from mule deer locations (xy.ppp) to the nearest river.

```{r eval=FALSE}
rivdist <- nncross(xy.ppp, segs2.psp)$dist

#Or identify segment number closest to each point
riv <- nearestsegment(xy.ppp,segs2.psp)
plot(segs2.psp, lwd=1)
plot(xy.ppp[1], add=TRUE, col="red")
plot(segs2.psp[riv[1]], add=TRUE, lwd=5)
plot(xy.ppp[10], add=TRUE, col="blue")
plot(segs2.psp[riv[10]], add=TRUE, lwd=5, col="blue")
```

18\. We can then summarize the distances in some meaningful way for analysis. Instead of representing distance to road as individual numerical values we can bin the distances in some categories we determine appropriate for our research objective.

```{r eval=FALSE}
br <- seq(0,1000,200)
lbl <- paste(head(br,-1),tail(br,-1),sep="-")
road.tbl <- table(cut(roaddist,breaks=br,labels=lbl))
Rdresults <- road.tbl/sum(road.tbl)
Rdresults

br1 <- seq(0,4000,500)
lbl1 <- paste(head(br1,-1),tail(br1,-1),sep="-")
river.tbl <- table(cut(rivdist,breaks=br1,labels=lbl1))
Rivresults <- river.tbl/sum(river.tbl)
Rivresults
```

19\. Or we can place each distance into a category or Bin for each deer

```{r eval=FALSE}
BinRoad <- bin(roaddist, nbins=5, method='content', labels=c('1','2','3','4','5'))
BinRoad2 <- cut(roaddist, 5, method='intervals', include.lowest=TRUE, labels=c('1','2','3'
  ,'4','5'))
table(BinRoad)

BinRivers <- bin(rivdist, nbins=5, method='content', labels=c('1','2','3','4','5'))
BinRivers <- cut(rivdist, 5, method='intervals', include.lowest=TRUE, labels=c('1','2','3'
    ,'4','5'))
table(BinRivers)

#Now use cbind function to add binned distances to muleys dataset.
#First I had to find a function to easily convert sf dataframe with geometry to dataframe that includes all data and columns for x and y:
#https://github.com/r-spatial/sf/issues/231
sfc_as_cols <- function(x, geometry, names = c("x","y")) {
  if (missing(geometry)) {
    geometry <- sf::st_geometry(x)
  } else {
    geometry <- rlang::eval_tidy(enquo(geometry), x)
  }
  stopifnot(inherits(x,"sf") && inherits(geometry,"sfc_POINT"))
  ret <- sf::st_coordinates(geometry)
  ret <- tibble::as_tibble(ret)
  stopifnot(length(names) == ncol(ret))
  x <- x[ , !names(x) %in% names]
  ret <- setNames(ret,names)
  dplyr::bind_cols(x,ret)
}

muleys2 <- sfc_as_cols(deer.albers, st_centroid(geometry))

Dist <- cbind(BinRoad,BinRivers)
muleys <- cbind(muleys2, roaddist, rivdist,Dist)
```