WG_Brecon_Fine.Rmd

---
title: "CUUPeat: Brecon Beacons, Wales - Fine-scale bare peat mapping"
date: "15/09/2020"
author: 'JNCC'
licence: 'MIT Licence'
output:
  html_document:
    df_print: paged
    css: style.css
    includes:
      before_body: header.html
  pdf_document: default
always_allow_html: yes
---

```{r setup, include=FALSE,echo=F}

knitr::opts_chunk$set(echo = TRUE)
# Built with r 3.6.0
library(sf)
library(purrr)
library(tmap)
library(dplyr)
library(raster)

#package dependencies
## purrr_0.3.4       tmap_3.0          sf_0.9-2         
## dplyr_0.8.5       raster_3.0-7     fasterize_1.0.2 
```

<div class="mycontent">

CUU Peat Project - Welsh Government - Brecon Beacons, Wales
Fine-scale bare peat mapping processing


## Site location

```{r sitemap,message=F,warning=F}

#site boundaries
site <- sf::st_read('Breacon_Site_Boundaries.shp',quiet=T)
#S2 AOI
S2 <- sf::st_read('Proposed_S2_AOI.shp',quiet=T)
#APGB AOI
APGB <- sf::st_read('Proposed_APGB_AOI.shp',quiet=T)

#ploy
tmap::tmap_mode("view")

tmap::tm_shape(site) +  tmap::tm_borders(alpha=0.5) +  tmap::tm_fill("green") +
  tmap::tm_shape(S2) +  tmap::tm_borders(alpha=1,col='blue') +
  tmap::tm_shape(APGB) +  tmap::tm_borders(alpha=1,col='pink')


```


## Preparing the APGB

For the Breacon Beacons site, there were APGB imagery available from:

* 2017-05-26
* 2018-05-07

```{r apgb_latest, include=T,message=F,warning=F,fig.show = "hold", out.width = "50%", fig.align = "default",eval=F}
## latest available data
dirpath <- "WG_BB/"

#run function to resample nir band from image b, and combine with image a. the nir argument specifys which band is the nir band in image b and the start and end arguments are for iterating between many data files. 

AOIs <- list.files(paste0(dirpath,"25cm Aerial Photo/"))
purrr::map(AOIs,.f=function(shape){
  comboAPGB(a.path=paste0(dirpath,"25cm Aerial Photo/",shape,"/"),
          b.path=paste0(dirpath,"50cm Colour Infrared/",gsub(shape,pattern="1_RGB",replacement="2_CIR"),"/"),
          nir=1,
          out.path=paste0(dirpath),
          start=1,end=NULL)
})

## rename files with dates
date_lookup <- 'CIR_Date_Grid.shp'

APGBdate(imgfolder=paste0(dirpath,'Combined/Shape1/'),
         lookup=date_lookup)
APGBdate(imgfolder=paste0(dirpath,'Combined/Shape2/'),
         lookup=date_lookup)
APGBdate(imgfolder=paste0(dirpath,'Combined/Shape3/'),
         lookup=date_lookup)

```


## Creating fine scale maps

### Subsetting the APGB tiles to the peat soils/site boundary layer 

```{r peatsoils,eval=F }
# prepare soil mask from peat soils data
peatsoils <- sf::st_read('CUU_WG_Peat_Training_Data.shp',quiet=T)

highprobpeat <- peatsoils %>% dplyr::filter(Peat=="Highly probable")

dir.create('processed')
sf::st_write(highprobpeat,'processed/fine_peat_soils.shp')

#site boundaries

# bare peat function
MaskPeatSoils <- function(img_folder,mask, delete=T){
  imgs <- list.files(img_folder, pattern='.tif', full.names = T )
  mask_sf <- sf::st_read(mask) %>% sf::st_transform(27700)
  #iterate through images
  purrr::map(imgs, .f=function(img){
    img_r <- raster::brick(img)

    #check if overlaps with the mask
    img_poly <- sf::st_as_sf(as(raster::extent(img_r), 'SpatialPolygons')) %>% sf::st_set_crs(27700)
    int_poly <- suppressWarnings(sf::st_intersection(mask_sf,img_poly))
    # if no overlap remove tile
    if(nrow(int_poly)==0){
      file.remove(img)
      print(paste0(basename(img), " done."))
    } else {
      # if overlap then mask tile to peat soil
    r <-raster::raster(ext = raster::extent(img_r),crs = raster::crs(img_r),res = raster::res(img_r))
    mask_r <- fasterize::fasterize(mask_sf, r)
    img_masked <-raster::overlay(img_r,mask_r,fun = function(x, y) {
      x[is.na(y[])] <- NA
      return(x)
    })
    raster::writeRaster(img_masked,img, overwrite=T)
    print(paste0(basename(img), " done."))
    }

  })

}
# apply the bare peat function to the different sites to only select areas where there is peat soil
dir_fold <- 'PeatlandCondition/'
# site 1
MaskPeatSoils(paste0(dir_fold,"Data/APGB/WG_BB/Combined/Shape1"),
              paste0(dir_fold,"Data/Under_licence/OFFICIAL_SENSITIVE_WGPeatSoils/processed/fine_peat_soils.shp"))
# site 2
MaskPeatSoils(paste0(dir_fold,"Data/APGB/WG_BB/Combined/Shape2"),
              paste0(dir_fold,"Data/Under_licence/OFFICIAL_SENSITIVE_WGPeatSoils/processed/fine_peat_soils.shp"))
# site 3
MaskPeatSoils(paste0(dir_fold,"Data/APGB/WG_BB/Combined/Shape3"),
              paste0(dir_fold,"Data/Under_licence/OFFICIAL_SENSITIVE_WGPeatSoils/processed/fine_peat_soils.shp"))
  
```

Threshold levels were explored using the peatpal app on areas with clear bare peat. These thresholds were then run across the whole of the site. The outputs were then examined in QGIS. Where there were areas where the thresholding rules did not seem to fit well (which might have happened if tiles came from different years for example), these files were examined again in peatpal to create new theshold levels.

```{r threshold, eval=F}
# load the barethresh function and the runIndices function

## threshold based on SN8116
#x1=nir,x2=brightness,x3=gli,x4=NDVI,x5=NDWI,x6=rb,x7=rg
thresh_fun <- function(x1,x2,x3,x4,x5,x6,x7) {
  ifelse(x1 <119.74 & x2 >70 & x2<110 & x3<0.03 & x4<0.07 & x5<0.1 & x6>1.25 & x7>0.01 & x7<1.12 , 1, NA) 
}

##site 1
img_folder <-paste0(dir_fold,'Data/APGB/WG_BB/Combined/Shape1/')
# threhold based on SN8217
#x1=nir,x2=brightness,x3=gli,x4=NDWI
thresh_fun <- function(x1,x2,x3,x4) {
  ifelse(x1 <107.25 & x2 >70 & x2<110 & x3<0.02 & x4<0.03, 1,NA) 
}
#run thresholding rule# all bands
barethresh(Img.path=img_folder, 
           out.path=img_folder,
           spec.bands=4, #if spectral bands are included in the thresholding, the band number
           ind.name=c("Brightness", "GLI","NDVI","NDWI","RB","RG"), #name of indices, should be in alphabetical order and will be used in this order in the thresholding function
           c.fun=thresh_fun, #thresholding function
           nir=4, r=1,g=2,b=3, #band numbers within the imagery
           start=1)

# site 1, tile SN8118
#x1=nir,x2=brightness,x3=gli,x4=NDWI
img_folder <- paste0(dir_fold,'Data/APGB/WG_BB/Combined/Shape1_SN8118/')
thresh_fun <- function(x1,x2,x3,x4) {
  ifelse(x1>75.08 & x1<118.15 & x2 >81 & x2<104 & x3>-0.03 & x3<0.02 & x4>-0.02 & x4<0.09, 1, NA) 
}
#run thresholding rule# all bands
barethresh(Img.path=img_folder, 
           out.path=img_folder,
           spec.bands=4, #if spectral bands are included in the thresholding, the band number
           ind.name=c("Brightness", "GLI","NDWI"), #name of indices, should be in alphabetical order and will be used in this order in the thresholding function
           c.fun=thresh_fun, #thresholding function
           nir=4, r=1,g=2,b=3, #band numbers within the imagery
           start=1)

# site 1, tile SN8218
#x1=nir,x2=brightness,x3=gli,x4=NDWI
img_folder <- paste0(dir_fold,'Data/APGB/WG_BB/Combined/Shape1_SN8218/')
thresh_fun <- function(x1,x2,x3,x4) {
  ifelse(x1>75.08 & x1<118.15 & x2 >81 & x2<104 & x3>-0.03 & x3<0.02 & x4>-0.02 & x4<0.09, 1, NA) 
}
#run thresholding rule# all bands
barethresh(Img.path=img_folder, 
           out.path=img_folder,
           spec.bands=4, #if spectral bands are included in the thresholding, the band number
           ind.name=c("Brightness", "GLI","NDWI"), #name of indices, should be in alphabetical order and will be used in this order in the thresholding function
           c.fun=thresh_fun, #thresholding function
           nir=4, r=1,g=2,b=3, #band numbers within the imagery
           start=1)

#site 2
img_folder <-paste0(dir_fold,'Data/APGB/WG_BB/Combined/Shape2/')
#thresholding rules, some bands
#run thresholding rule
barethresh(Img.path=img_folder, 
           out.path=img_folder,
           spec.bands=4, #if spectral bands are included in the thresholding, the band number
           ind.name=c("Brightness", "GLI","NDWI"), #name of indices, should be in alphabetical order and will be used in this order in the thresholding function
           c.fun=thresh_fun, #thresholding function
           nir=4, r=1,g=2,b=3, #band numbers within the imagery
           start=1)
# threshold based on 2018-05-07_SO0119.tif
#x1=nir,x2=brightness,x3=gli,x4=NDVI,x5=NDWI
thresh_fun <- function(x1,x2,x3,x4,x5) {
  ifelse(x1 <79.85 & x2 >60 & x2<80 & x3<0.01 & x4<0 & x5<0.09, 1,NA) 
}
#run thresholding rule# all bands
barethresh(Img.path=img_folder, 
           out.path=img_folder,
           spec.bands=4, #if spectral bands are included in the thresholding, the band number
           ind.name=c("Brightness", "GLI","NDVI","NDWI"), #name of indices, should be in alphabetical order and will be used in this order in the thresholding function
           c.fun=thresh_fun, #thresholding function
           nir=4, r=1,g=2,b=3, #band numbers within the imagery
           start=1)
# threshold based on 4_6_2017-05-26_SO0017.tif
#x1=nir,x2=brightness,x3=gli,x4=NDVI,x5=NDWI
thresh_fun <- function(x1,x2,x3,x4,x5) {
  ifelse(x1 <109.5 & x2 >70 & x2<96 & x3<0.01 & x4<0 & x5<0.2, 1,NA) 
}
#run thresholding rule# all bands
barethresh(Img.path=img_folder, 
           out.path=img_folder,
           spec.bands=4, #if spectral bands are included in the thresholding, the band number
           ind.name=c("Brightness", "GLI","NDVI","NDWI"), #name of indices, should be in alphabetical order and will be used in this order in the thresholding function
           c.fun=thresh_fun, #thresholding function
           nir=4, r=1,g=2,b=3, #band numbers within the imagery
           start=1)

#site 3
img_folder <-paste0(dir_fold,'Data/APGB/WG_BB/Combined/Shape3/')
## threshold based on SN2118
#x1=nir,x2=brightness,x3=gli,x4=NDVI,x5=NDWI
thresh_fun <- function(x1,x2,x3,x4,x5) {
  ifelse(x1 <90.24 & x2 >70 & x2<105 & x3<0.02 & x4<0 & x5<0.12, 1,NA) 
}
#run thresholding rule# all bands
barethresh(Img.path=img_folder, 
           out.path=img_folder,
           spec.bands=4, #if spectral bands are included in the thresholding, the band number
           ind.name=c("Brightness", "GLI","NDVI","NDWI"), #name of indices, should be in alphabetical order and will be used in this order in the thresholding function
           c.fun=thresh_fun, #thresholding function
           nir=4, r=1,g=2,b=3, #band numbers within the imagery
           start=1)
```

The final bare peat files were moved into the folders: 
'J:/GISprojects/EOMonitoringApplications/CUU_WP6_Projects/PeatlandCondition/WG_BreconBeacons/Fine_scale_mapping/Site1/Bare'
'J:/GISprojects/EOMonitoringApplications/CUU_WP6_Projects/PeatlandCondition/WG_BreconBeacons/Fine_scale_mapping/Site2/Bare'
'J:/GISprojects/EOMonitoringApplications/CUU_WP6_Projects/PeatlandCondition/WG_BreconBeacons/Fine_scale_mapping/Site3/Bare'

All the tiles converted in a binary classification; 0 if there was no bare peat present and 1 if there was bare peat present. This was converted to the 10m pixel level as a percentage cover per 10m of bare peat.

```{r, eval=F}
# for each site, create layers of 10m percentage peat cover #
### SITE 1
dir_fold <- 'PeatlandCondition/'

dirpath <- paste0(dir_fold,'Data/APGB/WG_BB/Combined/Shape1/')
#convert into percentage covers
pcov_classify(img_path=paste0(dirpath,'masked/'), 
              bare_path=paste0(dirpath,'bare/'),
              out_path=dirpath,
              polymask='fine_peat_soils.shp')
### SITE 2
dirpath <- paste0(dir_fold,'Data/APGB/WG_BB/Combined/Shape2/')
#convert into percentage covers
pcov_classify(img_path=paste0(dirpath,'masked/'), 
              bare_path=paste0(dirpath,'bare/'),
              out_path=dirpath,
              polymask='fine_peat_soils.shp')
### SITE 3
dirpath <- paste0(dir_fold,'Data/APGB/WG_BB/Combined/Shape3/')
#convert into percentage covers
pcov_classify(img_path=paste0(dirpath,'masked/'), 
              bare_path=paste0(dirpath,'bare/'),
              out_path=dirpath,
              polymask='fine_peat_soils.shp')


# mosaic tiles
dirpath <- paste0(dir_fold,'Data/APGB/WG_BB/Combined/')
### 2017
pcov <- list.files(paste0(dirpath,'perc_cov2017/'),full.names=T)
pcov_list <- purrr::map(pcov,raster)
names(pcov_list)<-NULL
pcov_list$fun <- mean
pcov_mosaic <- do.call(raster::mosaic,c(pcov_list,progress="window"))
raster::writeRaster(pcov_mosaic,paste0(dirpath,'bare_pcov_mosaic_2017.tif'))
### 2018
pcov <- list.files(paste0(dirpath,'perc_cov2018/'),full.names=T)
pcov_list <- purrr::map(pcov,raster)
names(pcov_list)<-NULL
pcov_list$fun <- mean
pcov_mosaic <- do.call(raster::mosaic,c(pcov_list,progress="window"))
raster::writeRaster(pcov_mosaic,paste0(dirpath,'bare_pcov_mosaic_2018.tif'))


```