
---
title: "Updated_SurfInt"
author: "Stephanie Burgess"
date: "06/27/2023"
output:
  pdf_document: default
  html_document: default
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

##The following code loads shapefiles for Industry infrastructure (all files must be downloaded into the working directory set in R); subsets the shapefile by occupancy rate in tenths (i.e., 0.1, 0.2, 0.3 etc); iteratively measures the area in each occupancy rate bin; calculates the coastal v. inland area; and outputs the area to a csv.
##update Jan 2023 included edits to accomodate runs that have <100% occupancy throughout the project.
##update June 2023 included the migration from sp-based packages to sf-based packages.
##update Sept 2023 included added narratives. 

##This code was created by Stephanie Burgess, USFWS, Alaska Region, FES-MMM. If you have any questions regarding the applicability of this analysis method, please contact Dr. Burgess at stephanie_burgess@fws.gov. 

###Add the libraries

```{r}

library(sf)
library(ggplot2)
library(rmapshaper)
```


### Set occupancy rates

```{r}
one.occ=0.1
two.occ=0.2
three.occ=0.3
four.occ=0.4
five.occ=0.5
six.occ=0.6
seven.occ=0.7
eight.occ=0.8
nine.occ=0.9
always.occ=1
```


### Read in the spatial objects and project to Albers Equal Area Projection. The code below relies upon shapefiles. All.buff should be a shapefile that is the merged result of buffering all spatial features (i.e., roads, pipelines, pads) by 1610 meters, withOUT dissolving into one buffer. Further, the file should include occupancy rates for every month of the year. Occupancy fields should be titled "MonYY" (e.g., Aug24). Occupancy rates should represent the amount of time during the month that humans are present at the location (this is a presence/absence metric, and does not reflect the number of people at a site). For example, if humans will be present every day of the month for 12 hours/ day, the occupancy rate would be 0.5. Similarly, if humans will be present 15 days out of 30 for 24 hours/day day, the occupancy rate would also be 0.5. Operators may choose to be as detailed with their occupancy estimates as they would like, however please note after averaging the occupancy for the season rates will be rounded up to the nearest 0.1 for computational tractability. Coastal and inland ITR areas were created by FWS and can be supplied as needed. 

```{r}

all.buff<-read_sf("PATH TO YOUR FILES/surface_analysis_input.shp")
all.buff<-st_transform(all.buff, st_crs(3338))
coast<-read_sf("PATH TO YOUR FILES/blm25_coast.shp")
coast<-st_transform(coast, st_crs(3338))
inland<-read_sf("PATH TO YOUR FILES/blm25_inland.shp")
inland<-st_transform(inland, st_crs(3338))
```


### Plot to make sure it looks right

```{r}
ggplot()+
     geom_sf(data=inland)+
     geom_sf(data=coast)+
     geom_sf(data=all.buff)+
     coord_sf(crs=st_crs(3338))
```

## Pick the season and year you want to analyze, and create seasonal occupancy rates. Please note the open season is 5 months, thus the denominator should always be 5, and the ice season is 7 months with a denominator of 7. 

```{r}
#year.seas<-all.buff$ice25
#label<-"ice25"

#all.buff$Open25<-((all.buff$Jul25)/5)
all.buff$ice25<-((all.buff$Jan25+all.buff$Feb25+all.buff$March25+all.buff$April25+all.buff$May25)/7)


year.seas<-all.buff$ice25
label<-"ice25"
project_name<-"blm_25"
```

### Subset the shapefile by occupancy rates for each season/year
```{r}
all.one<-subset(all.buff, year.seas>0 & year.seas<=0.1)
all.two<-subset(all.buff, year.seas>0.1 & year.seas<=0.2)
all.three<-subset(all.buff, year.seas>0.2 & year.seas<=0.3)
all.four<-subset(all.buff, year.seas>0.3 & year.seas<=0.4)
all.five<-subset(all.buff, year.seas>0.4 & year.seas<=0.5)
all.six<-subset(all.buff, year.seas>0.5 & year.seas<=0.6)
all.seven<-subset(all.buff, year.seas>0.6 & year.seas<=0.7)
all.eight<-subset(all.buff, year.seas>0.7 & year.seas<=0.8)
all.nine<-subset(all.buff, year.seas>0.8 & year.seas<=0.9)
all.always<-subset(all.buff, year.seas>0.9 & year.seas<=1)
```

### We now need to subtract the areas from one another by occupancy rate. These are the series of files that will be subtracted:

```{r}
from.one<-subset(all.buff, year.seas>0.1)
from.two<-subset(all.buff, year.seas>0.2)
from.three<-subset(all.buff, year.seas>0.3)
from.four<-subset(all.buff, year.seas>0.4)
from.five<-subset(all.buff, year.seas>0.5)
from.six<-subset(all.buff, year.seas>0.6)
from.seven<-subset(all.buff, year.seas>0.7)
from.eight<-subset(all.buff, year.seas>0.8)
from.nine<-subset(all.buff, year.seas>0.9)
```


### Get the difference between each file and the ones subtracted from it

```{r}

if((nrow(from.one)>0) && (nrow(all.one)>0)){just.one=ms_erase(all.one,from.one)}#if "from.one"is NOT null (ie there are areas that have greater occupancy than one to take away from it)then this line will return all the one bin with areas with greater occupancy taken from it (because they have already been counted)
if((nrow(from.one)==0)||(nrow(all.one)==0)){just.one=all.one}#if there is nothing with a greater occupancy than one to take away, then it will just return all.one

if((nrow(from.two)>0) && (nrow(all.two)>0)){just.two=ms_erase(all.two,from.two)}
if((nrow(from.two)==0)||(nrow(all.two)==0)){just.two=all.two}

if((nrow(from.three)>0) && (nrow(all.three)>0)){just.three=ms_erase(all.three,from.three)}
if((nrow(from.three)==0)||(nrow(all.three)==0)){just.three=all.three}

if((nrow(from.four)>0) && (nrow(all.four)>0)){just.four=ms_erase(all.four,from.four)}
if((nrow(from.four)==0)||(nrow(all.four)==0)){just.four=all.four}

if((nrow(from.five)>0) && (nrow(all.five)>0)){just.five=ms_erase(all.five,from.five)}
if((nrow(from.five)==0)||(nrow(all.five)==0)){just.five=all.five}

if((nrow(from.six)>0) && (nrow(all.six)>0)){just.six=ms_erase(all.six,from.six)}
if((nrow(from.six)==0)||(nrow(all.six)==0)){just.six=all.six}

if((nrow(from.seven)>0) && (nrow(all.seven)>0)){just.seven=ms_erase(all.seven,from.seven)}
if((nrow(from.seven)==0)||(nrow(all.seven)==0)){just.seven=all.seven}

if((nrow(from.eight)>0) && (nrow(all.eight)>0)){just.eight=ms_erase(all.eight,from.eight)}
if((nrow(from.eight)==0)||(nrow(all.eight)==0)){just.eight=all.eight}

if((nrow(from.nine)>0) && (nrow(all.nine)>0)){just.nine=ms_erase(all.nine,from.nine)}
if((nrow(from.nine)==0)||(nrow(all.nine)==0)){just.nine=all.nine}
```


### Compute coastal and inland areas

```{r}
test.all.c<-st_intersects(st_union(all.always),st_union(coast))
test.all.i<-st_intersects(st_union(all.always),st_union(inland))
if(nrow(all.always)>0){
if(sum(unlist(test.all.c))>0) {
  always.area.coast=as.numeric(sum(st_area(st_intersection(st_union(all.always),st_union(coast)))))/1000^2}
if(sum(unlist(test.all.i))>0) {
always.area.inland=as.numeric(sum(st_area(st_intersection(st_union(all.always),st_union(inland)))))/1000^2}}

test.nine.c<-st_intersects(st_union(just.nine),st_union(coast))
test.nine.i<-st_intersects(st_union(just.nine),st_union(inland))
if(nrow(just.nine)>0){
if(sum(unlist(test.nine.c))>0) {
  nine.area.coast=as.numeric(sum(st_area(st_intersection(st_union(just.nine),st_union(coast)))))/1000^2}
if(sum(unlist(test.nine.i))>0) {
nine.area.inland=as.numeric(sum(st_area(st_intersection(st_union(just.nine),st_union(inland)))))/1000^2}}

test.eight.c<-st_intersects(st_union(just.eight),st_union(coast))
test.eight.i<-st_intersects(st_union(just.eight),st_union(inland))
if(nrow(just.eight)>0){
if(sum(unlist(test.eight.c))>0) {
  eight.area.coast=as.numeric(sum(st_area(st_intersection(st_union(just.eight),st_union(coast)))))/1000^2}
if(sum(unlist(test.eight.i))>0) {
eight.area.inland=as.numeric(sum(st_area(st_intersection(st_union(just.eight),st_union(inland)))))/1000^2}}

test.seven.c<-st_intersects(st_union(just.seven),st_union(coast))
test.seven.i<-st_intersects(st_union(just.seven),st_union(inland))
if(nrow(just.seven)>0){
if(sum(unlist(test.seven.c))>0) {
  seven.area.coast=as.numeric(sum(st_area(st_intersection(st_union(just.seven),st_union(coast)))))/1000^2}
if(sum(unlist(test.seven.i))>0) {
seven.area.inland=as.numeric(sum(st_area(st_intersection(st_union(just.seven),st_union(inland)))))/1000^2}}

test.six.c<-st_intersects(st_union(just.six),st_union(coast))
test.six.i<-st_intersects(st_union(just.six),st_union(inland))
if(nrow(just.six)>0){
if(sum(unlist(test.six.c))>0) {
  six.area.coast=as.numeric(sum(st_area(st_intersection(st_union(just.six),st_union(coast)))))/1000^2}
if(sum(unlist(test.six.i))>0) {
six.area.inland=as.numeric(sum(st_area(st_intersection(st_union(just.six),st_union(inland)))))/1000^2}}

test.five.c<-st_intersects(st_union(just.five),st_union(coast))
test.five.i<-st_intersects(st_union(just.five),st_union(inland))
if(nrow(just.five)>0){
if(sum(unlist(test.five.c))>0) {
  five.area.coast=as.numeric(sum(st_area(st_intersection(st_union(just.five),st_union(coast)))))/1000^2}
if(sum(unlist(test.five.i))>0) {
five.area.inland=as.numeric(sum(st_area(st_intersection(st_union(just.five),st_union(inland)))))/1000^2}}

test.four.c<-st_intersects(st_union(just.four),st_union(coast))
test.four.i<-st_intersects(st_union(just.four),st_union(inland))
if(nrow(just.four)>0){
if(sum(unlist(test.four.c))>0) {
  four.area.coast=as.numeric(sum(st_area(st_intersection(st_union(just.four),st_union(coast)))))/1000^2}
if(sum(unlist(test.four.i))>0) {
four.area.inland=as.numeric(sum(st_area(st_intersection(st_union(just.four),st_union(inland)))))/1000^2}}

test.three.c<-st_intersects(st_union(just.three),st_union(coast))
test.three.i<-st_intersects(st_union(just.three),st_union(inland))
if(nrow(just.three)>0){
if(sum(unlist(test.three.c))>0) {
  three.area.coast=as.numeric(sum(st_area(st_intersection(st_union(just.three),st_union(coast)))))/1000^2}
if(sum(unlist(test.three.i))>0) {
three.area.inland=as.numeric(sum(st_area(st_intersection(st_union(just.three),st_union(inland)))))/1000^2}}

test.two.c<-st_intersects(st_union(just.two),st_union(coast))
test.two.i<-st_intersects(st_union(just.two),st_union(inland))
if(nrow(just.two)>0){
if(sum(unlist(test.two.c))>0) {
  two.area.coast=as.numeric(sum(st_area(st_intersection(st_union(just.two),st_union(coast)))))/1000^2}
if(sum(unlist(test.two.i))>0) {
two.area.inland=as.numeric(sum(st_area(st_intersection(st_union(just.two),st_union(inland)))))/1000^2}}

test.one.c<-st_intersects(st_union(just.one),st_union(coast))
test.one.i<-st_intersects(st_union(just.one),st_union(inland))
if(nrow(just.one)>0){
if(sum(unlist(test.one.c))>0) {
  one.area.coast=as.numeric(sum(st_area(st_intersection(st_union(just.one),st_union(coast)))))/1000^2}
if(sum(unlist(test.one.i))>0) {
one.area.inland=as.numeric(sum(st_area(st_intersection(st_union(just.one),st_union(inland)))))/1000^2}}

```

### Check that everything adds up. This may deliver FALSE if the values are even slightly different. Inspect the areas if this occurs, they may be extremely close. In that event, use a visual inspection of the shapefiles to identify the issue, or you may make the decision to use the higher value conservatively. 

```{r}
occupied<-st_union(subset(all.buff, year.seas>0))
tot1=st_intersection(occupied,st_union(coast))
round(((ifelse(exists("one.area.coast"),one.area.coast,0))
+(ifelse(exists("two.area.coast"),two.area.coast,0))              
+(ifelse(exists("three.area.coast"),three.area.coast,0))
+(ifelse(exists("four.area.coast"),four.area.coast,0))             
+(ifelse(exists("five.area.coast"),five.area.coast,0))
+(ifelse(exists("six.area.coast"),six.area.coast,0))
+(ifelse(exists("seven.area.coast"),seven.area.coast,0))
+(ifelse(exists("eight.area.coast"),eight.area.coast,0))
+(ifelse(exists("nine.area.coast"),nine.area.coast,0))
+(ifelse(exists("always.area.coast"),always.area.coast,0))),1)==round((as.numeric(sum(st_area(tot1)))/1000^2),1) #round due to numerical error in functions

tot2=st_intersection(occupied,st_union(inland))
round(((ifelse(exists("one.area.inland"),one.area.inland,0))
+(ifelse(exists("two.area.inland"),two.area.inland,0))              
+(ifelse(exists("three.area.inland"),three.area.inland,0))
+(ifelse(exists("four.area.inland"),four.area.inland,0))             
+(ifelse(exists("five.area.inland"),five.area.inland,0))
+(ifelse(exists("six.area.inland"),six.area.inland,0))
+(ifelse(exists("seven.area.inland"),seven.area.inland,0))
+(ifelse(exists("eight.area.inland"),eight.area.inland,0))
+(ifelse(exists("nine.area.inland"),nine.area.inland,0))
+(ifelse(exists("always.area.inland"),always.area.inland,0))),1)==round((as.numeric(sum(st_area(tot2)))/1000^2),1) #round due to numerical error in functions

ggplot()+
     geom_sf(data=inland, color="blue")+
     geom_sf(data=coast, color="red")+
     geom_sf(data=tot1, fill="purple")+
     geom_sf(data=tot2, fill="yellow")+
     coord_sf(crs=st_crs(3338))

```


### Compile and write output

```{r}
(area.occ.df=data.frame(area=c((ifelse(exists("one.area.coast"),one.area.coast,0)),
                               (ifelse(exists("two.area.coast"),two.area.coast,0)), (ifelse(exists("three.area.coast"),three.area.coast,0)),
                               (ifelse(exists("four.area.coast"),four.area.coast,0)),             
                               (ifelse(exists("five.area.coast"),five.area.coast,0)),
                               (ifelse(exists("six.area.coast"),six.area.coast,0)),
                               (ifelse(exists("seven.area.coast"),seven.area.coast,0)),
                               (ifelse(exists("eight.area.coast"),eight.area.coast,0)),
                               (ifelse(exists("nine.area.coast"),nine.area.coast,0)),
                               (ifelse(exists("always.area.coast"),always.area.coast,0)),
                               (ifelse(exists("one.area.inland"),one.area.inland,0)),
                               (ifelse(exists("two.area.inland"),two.area.inland,0)), (ifelse(exists("three.area.inland"),three.area.inland,0)),
                               (ifelse(exists("four.area.inland"),four.area.inland,0)),            
                               (ifelse(exists("five.area.inland"),five.area.inland,0)),
                               (ifelse(exists("six.area.inland"),six.area.inland,0)),
                               (ifelse(exists("seven.area.inland"),seven.area.inland,0)),
                               (ifelse(exists("eight.area.inland"),eight.area.inland,0)),
                               (ifelse(exists("nine.area.inland"),nine.area.inland,0)),
                               (ifelse(exists("always.area.inland"),always.area.inland,0)),
                        occ=c(one.occ, two.occ, three.occ, four.occ, five.occ, six.occ, seven.occ, eight.occ, nine.occ, always.occ, one.occ, two.occ, three.occ, four.occ, five.occ, six.occ, seven.occ, eight.occ, nine.occ, always.occ),
                        region=c(rep('coast',10),
                                 rep('inland',10)))))
write.csv(area.occ.df,paste0(project_name,label,'.csv'))


```          

### **IMPORTANT:** Remove all of the objects we just created so they are not accidently used in subsequent runs. 

```{r}
rm(list=ls())

```
