### Modeling range of Varecia, the Black-and-White Lemur (with Toni Lyn Morelli)

### Adam B. Smith | Missouri Botanical Garden | adamDOTsmithATmobotDOTorg | 2018-02

### This code is version 6 of attempts to model the exposure of two lemurs, Varecia variegata and Varecia rubra, to anticipated climate change and deforestation in Madagascar. The analysis depends on two models, one of deforestation and an ecological niche model using forest cover and climate as predictors.

# source('C:/Ecology/Drive/Research/Varecia (Toni Lyn Morelli)/Versions 06/Code/Varecia Ver 06 Forest Modeling.r')

# source('E:/Ecology/Drive/Research/Varecia (Toni Lyn Morelli)/Versions 06/Code/Varecia Ver 06 Forest Modeling.r')

# source('H:/Global Change Program/Research/Varecia (Toni Lyn Morelli)/Versions 06/Code/Varecia Ver 06 Forest Modeling.r')

setwd('C:/ecology/Drive/Research/Varecia (Toni Lyn Morelli)/Versions 06')

# setwd('E:/ecology/Drive/Research/Varecia (Toni Lyn Morelli)/Versions 06')

# setwd('H:/Global Change Program/Research/Varecia (Toni Lyn Morelli)/Versions 06')

################

### CONTENTS ###

################

### SETUP ###

### libraries, functions, and definitions ###

### define study region: humid forest ###

### crop forest maps to study region ###

### MODELING DEFORESTATION ###

### collate data for modeling deforestation ###

### estimate deforestation VOLUME (vs rate) ###

### Vielledent model of deforestation RATE ###

### Vielledent model of deforestation LOCATION ###

### evaluate deforestation LOCATION models ###

### project deforestation into future ###

### create maps of future forest cover assuming no additional loss in protected areas after 2014 ###

### calculate future forest fragmentation class ###

### compile forest cover and fragmentation statistics ###

### create display maps of forest cover ###

### compile GIS files for faster plotting ###

### create display maps of just Madagascar for presentations ###

### create display maps of PAST forest cover for presentations ###

### create display maps of FUTURE forest cover for presentations ###

### create side-by-side display maps of FUTURE forest cover for animations ###

### create display maps of forest cover in MOBOT reserves for presentations ###

### create 3D display maps of forest cover ###

### create display maps of deforestation probability ###

### create display maps of forest fragmentation class ###

#############################################

### libraries, functions, and definitions ###

#############################################

memory.limit(memory.limit() * 2^30)

rm(list=ls())

gc()

library(omnibus)

library(enmSdm)

library(legendary)

library(fasterRaster)

library(sp)

library(graticule)

library(raster)

library(parallel)

library(dismo)

library(rgeos)

# library(phcfM)

library(purrr)

library(doBy)

# library(rayshader)

library(rgl)

library(scales)

library(tictoc)

say(date())

###############

### options ###

###############

options(stringsAsFactors=FALSE)

rasterOptions(format='GTiff', overwrite=TRUE, tmpdir='C:/ecology/!Scratch/_raster')

#################

### variables ###

#################

madEaProj <- '+init=epsg:32738'

bioclims <- c(4, 10, 15, 16) # BIOCLIM variables to use

gcms <- c('CCSM4', 'MIROC-ESM-CHEM', 'GISS-E2-R') # GCMs to use

longLat <- c('LONG_WGS84', 'LAT_WGS84') # BIOCLIM variables to use

neighs <- c(3, 5) # neighborhood size (in # of cells) to use to calculate forest loss

grassDir <- c('C:/OSGeo4W64/', 'grass-7.4.1', 'osgeo4W')

### extents of insets for plots of maps

madFocus1 <- extent(929471, 1043462, 8234420, 8358325)

madFocus2 <- extent(850269, 944754, 7949353, 8048591)

madFocus3 <- extent(687756, 751211.0, 7433385, 7553616)

madFocus1 <- as(madFocus1, 'SpatialPolygons')

madFocus2 <- as(madFocus2, 'SpatialPolygons')

madFocus3 <- as(madFocus3, 'SpatialPolygons')

projection(madFocus1) <- projection(madEaProj)

projection(madFocus2) <- projection(madEaProj)

projection(madFocus3) <- projection(madEaProj)

#################

### functions ###

#################

# equation for annualized forest loss rate from Puyravaud, J-P. 2003 Standardizing the calculation of the annual rate of deforestation. Forest Ecology and Management 177:593-596.

# t1, t1 Time 1 and 2

# A1, A2 Area of forest in times 1 and 2

compoundInterestLaw <- function(t1, t2, A1, A2) (1 / (t2 - t1)) * log(A2 / A1)

### convert singlepart spatial polygons to multipart

multiPart <- function(shape) {

}

# say('#########################################')

# say('### define study region: humid forest ###')

# say('#########################################')

# say('Ecoregions of Madagascar from citations in Vieilledent.')

# say('### UTM 38S ###')

# # get humid forest ecoregion

# ecoregions_utm38s <- shapefile('./Data/Ecoregions of Madagascar - Citations in Vieilledent et al 2018 Fig 1/madagascar_ecoregion_tenaizy_38s')

# ecoregions_utm38s <- gBuffer(ecoregions_utm38s, width=0, byid=TRUE) # fixed problem with self-intersection in shapefile

# humidForest_utm38s <- ecoregions_utm38s[ecoregions_utm38s$nature == 'Ecoregion Est', ]

# humidForest_utm38s <- sp::spTransform(humidForest_utm38s, CRS(madEaProj))

# # remove holes

# p4s <- projection(humidForest_utm38s)

# humidNoHole <- SpatialPolygons(list(Polygons(list(humidForest_utm38s@polygons[[1]]@Polygons[[1]]), ID=1)))

# humidForest_utm38s <- SpatialPolygonsDataFrame(humidNoHole, data=humidForest_utm38s@data, match.ID=FALSE)

# projection(humidForest_utm38s) <- p4s

# # remove disjunct areas in north

# humidForest_utm38s <- multiPart(humidForest_utm38s)

# humidForest_utm38s$area_km2 <- gArea(humidForest_utm38s, byid=TRUE) / (1000^2)

# humidForest_utm38s <- humidForest_utm38s[which.max(humidForest_utm38s$area_km2), ]

# # buffer to add all presences/absences for Varecia

# varecia <- read.csv('./Data/Varecia/02 Master__Varecia_Point_Data_2017 - Removed and Corrected Badly Georeferenced Sites and Zahamena site by EEL on 3-10-2000.csv')

# varecia <- SpatialPointsDataFrame(varecia[ , longLat], data=varecia, proj4string=getCRS('wgs84', TRUE))

# # get distance between outlying site and ecoregion

# varecia_utm38s <- sp::spTransform(varecia, CRS(madEaProj))

# distToHumid <- gDistance(varecia_utm38s, humidForest_utm38s, byid=TRUE)

# maxDist <- max(distToHumid)

# # add buffer to humid forest ecoregion and crop

# bufferSize <- 3 * maxDist

# humidForestBuffer_utm38s <- gBuffer(humidForest_utm38s, width=bufferSize)

# humidForestBuffer_utm38s <- crop(humidForestBuffer_utm38s, ecoregions_utm38s)

# # create mask

# forest <- raster('./Data/Forest - Vieilledent et al 2018/ORIGINALS/for2000.tif')

# proj4string(forest) <- madEaProj

# forest <- crop(forest, humidForestBuffer_utm38s)

# humidForestBufferMask_utm38s <- fasterRasterize(humidForestBuffer_utm38s, forest, grassDir=grassDir)

# names(humidForestBufferMask_utm38s) <- 'humidForestBufferMask_utm38s'

# # country polygon

# madagascar_wgs84 <- getData('GADM', country='MDG', level=0, path='./Study Region & Masks')

# madagascar_wgs84 <- gadm[gadm$NAME_0 == 'Madagascar', ]

# madagascar_utm38s <- sp::spTransform(madagascar_wgs84, CRS(madEaProj))

# # Jesus saves

# dirCreate('./Study Region & Masks/UTM 38S 30-m Resolution')

# save(humidForest_utm38s, file='./Study Region & Masks/UTM 38S 30-m Resolution/Eastern Humid Forest Polygon.RData')

# save(humidForestBuffer_utm38s, file='./Study Region & Masks/UTM 38S 30-m Resolution/Eastern Humid Forest Polygon Buffer.RData')

# writeRaster(humidForestBufferMask_utm38s, './Study Region & Masks/UTM 38S 30-m Resolution/humidForestBufferMask_utm38s', datatype='INT1U')

# save(madagascar_utm38s, file='./Study Region & Masks/UTM 38S 30-m Resolution/Madagascar from GADM.RData')

# say('### WGS 84 30 arcsec resolution ###')

# # project to WGS84

# humidForest_wgs84 <- sp::spTransform(humidForest_utm38s, getCRS('wgs84', TRUE))

# humidForestBuffer_wgs84 <- sp::spTransform(humidForestBuffer_utm38s, getCRS('wgs84', TRUE))

# # create mask

# wc <- raster('D:/Ecology/Climate/WORLDCLIM Ver 2 Rel June 1 2016/30 arcsec 1970 to 2000/wc2.0_bio_30s_01.tif')

# wc <- crop(wc, humidForestBuffer_wgs84)

# humidForestBufferMask_wgs84 <- rasterize(humidForestBuffer_wgs84, wc)

# humidForestBufferMask_wgs84 <- 1 + 0 * humidForestBufferMask_wgs84

# names(humidForestBufferMask_wgs84) <- 'humidForestBufferMask_wgs84'

# # Jesus saves

# dirCreate('./Study Region & Masks/WGS84 30-arcsec Resolution')

# save(humidForest_wgs84, file='./Study Region & Masks/WGS84 30-arcsec Resolution/Eastern Humid Forest Polygon.RData')

# save(humidForestBuffer_wgs84, file='./Study Region & Masks/WGS84 30-arcsec Resolution/Eastern Humid Forest Polygon Buffer.RData')

# writeRaster(humidForestBufferMask_wgs84, './Study Region & Masks/WGS84 30-arcsec Resolution/humidForestBufferMask_wgs84', datatype='INT1U')

# save(madagascar_wgs84, file='./Study Region & Masks/WGS84 30-arcsec Resolution/Madagascar from GADM.RData')

# say('########################################')

# say('### crop forest maps to study region ###')

# say('########################################')

# years <- c(2000, 2005, 2010, 2014)

# humidForestBufferMask_utm38s <- raster('./Study Region & Masks/UTM 38S 30-m Resolution/humidForestBufferMask_utm38s.tif')

# for (year in years) {

# say(year)

# rast <- raster(paste0('./Data/Forest - Vieilledent et al 2018/ORIGINALS/for', year, '.tif'))

# rast <- crop(rast, humidForestBufferMask_utm38s)

# rast <- rast * humidForestBufferMask_utm38s

# names(rast) <- paste0('forest', year)

# writeRaster(rast, paste0('./Data/Forest - Vieilledent et al 2018/forest', year), datatype='INT1U')

# rm(rast); gc()

# }

say('###############################################')

say('### collate data for modeling deforestation ###')

say('###############################################')

humidForestBufferMask_utm38s <- raster('./Study Region & Masks/UTM 38S 30-m Resolution/humidForestBufferMask_utm38s.tif')

say('### long/lat rasters')

ll <- longLatRasters(humidForestBufferMask_utm38s)

ll <- ll * humidForestBufferMask_utm38s

names(ll) <- c('longitude', 'latitude')

dirCreate('./Data/Longitude & Latitude Rasters')

writeRaster(ll, './Data/Longitude & Latitude Rasters/longLat_utm38s', datatype='INT4S')

say('### distance to nearest settlement')

places <- read.csv('./Data/Geographic Places (NIMA from DIVA-GIS)/00 Places.csv')

populatedPlaces <- places[places$F_DESIG == 'PPL', ]

populatedPlaces <- SpatialPointsDataFrame(cbind(populatedPlaces$LONG, populatedPlaces$LAT), data=populatedPlaces, proj4string=getCRS('wgs84', TRUE))

populatedPlaces <- sp::spTransform(populatedPlaces, CRS(madEaProj))

rast <- fasterVectToRastDistance(rast=humidForestBufferMask_utm38s, vect=populatedPlaces, grassDir=grassDir)

rast <- round(rast)

rast <- humidForestBufferMask_utm38s * rast

names(rast) <- 'distToSettlementMeters_utm38s'

writeRaster(rast, './Data/Geographic Places (NIMA from DIVA-GIS)/distToSettlementMeters_utm38s', datatype='INT4S')

say('### distance to major roads')

vect <- shapefile('./Data/Roads/MDG_roads')

vect <- sp::spTransform(vect, CRS(madEaProj))

out <- fasterVectToRastDistance(vect=vect, rast=humidForestBufferMask_utm38s, metric='euclidean', meters=TRUE, grassDir=grassDir)

out <- out * humidForestBufferMask_utm38s

out <- round(out)

names(out) <- 'distToNearestRoadMeters_utm38s'

writeRaster(out, './Data/Roads/distToNearestRoadMeters_utm38s', datatype='INT4S')

say('### distance to any water body (marine/inland)')

### lakes

say('lakes')

vect <- shapefile('./Data/Water Bodies/MDG_water_areas_dcw')

vect <- sp::spTransform(vect, CRS(projection(madEaProj)))

vect <- crop(vect, humidForestBufferMask_utm38s)

out <- fasterVectToRastDistance(rast=humidForestBufferMask_utm38s, vect=vect, metric='euclidean', meters=TRUE, grassDir=grassDir)

out <- out * humidForestBufferMask_utm38s

out <- round(out)

names(out) <- 'distToNearestLakeMeters_utm38s'

writeRaster(out, './Data/Water Bodies/distToNearestLakeMeters_utm38s', datatype='INT4S')

### rivers

say('rivers')

vect <- shapefile('./Data/Water Bodies/MDG_water_lines_dcw')

vect <- sp::spTransform(vect, CRS(projection(madEaProj)))

vect <- crop(vect, humidForestBufferMask_utm38s)

out <- fasterVectToRastDistance(rast=humidForestBufferMask_utm38s, vect=vect, metric='euclidean', meters=TRUE, grassDir=grassDir)

out <- out * humidForestBufferMask_utm38s

out <- round(out)

names(out) <- 'distToNearestRiverMeters_utm38s'

writeRaster(out, './Data/Water Bodies/distToNearestRiverMeters_utm38s', datatype='INT4S')

say('### distance to any inland water ###')

water <- raster::stack(c(

'./Data/Water Bodies/distToNearestLakeMeters_utm38s.tif',

'./Data/Water Bodies/distToNearestRiverMeters_utm38s.tif'

))

beginCluster(7)

out <- calc(water, min)

endCluster()

names(out) <- 'distToNearestInlandWaterMeters_utm38s'

writeRaster(out, './Data/Water Bodies/distToNearestInlandWaterMeters_utm38s', datatype='INT4S')

say('### distance to coast')

load('./Study Region & Masks/UTM 38S 30-m Resolution/Madagascar from GADM.RData')

madMask_utm38s <- fasterRasterize(vect=madagascar_utm38s, rast=humidForestBufferMask_utm38s, grassDir=grassDir)

# calculate distance to nearest coast

out <- fasterRastDistance(madMask_utm38s, metric='euclidean', meters=TRUE, fillNAs=FALSE, grassDir=grassDir)

out <- out * humidForestBufferMask_utm38s

out <- round(out)

names(out) <- 'distToNearestCoastMeters_utm38s'

writeRaster(out, './Data/Water Bodies/distToNearestCoastMeters_utm38s', datatype='INT4S')

say('### distance to any water')

distToWaterBody <- raster::stack(c(

'./Data/Water Bodies/distToNearestLakeMeters_utm38s.tif',

'./Data/Water Bodies/distToNearestRiverMeters_utm38s.tif',

'./Data/Water Bodies/distToNearestCoastMeters_utm38s.tif')

)

out <- min(distToWaterBody)

names(out) <- 'distToNearestWaterMeters_utm38s'

writeRaster(out, './Data/Water Bodies/distToNearestWaterBodyMeters_utm38s')

say('### topography from GMTED2010')

### elevation

elev <- fasterProjectRaster('D:/Ecology/Topography/GMTED2010/7pt5_arcsec/30S030E_20101117_gmted_mea075.tif', template=humidForestBufferMask_utm38s, method='bilinear', grassDir=grassDir)

### slope

slope <- fasterTerrain(elev, slope=TRUE, slopeUnits='degrees', grassDir=grassDir)

slope <- slope * humidForestBufferMask_utm38s

slope <- round(slope)

names(slope) <- 'slope'

writeRaster(slope, './Data/Topography - GMTED2010/slopeGmted2010_utm38s', dataType='INT1U')

### elevation (crop & save)

elev <- elev * humidForestBufferMask_utm38s

elev <- round(elev)

names(elev) <- 'elevation_gmted2010'

dirCreate('./Data/Topography - GMTED2010')

writeRaster(elev, './Data/Topography - GMTED2010/elevationGmted2010_utm38s', dataType='INT2S')

say('### forest layers ###')

forest2000 <- raster('./Data/Forest - Vieilledent et al 2018/forest2000.tif')

forest2005 <- raster('./Data/Forest - Vieilledent et al 2018/forest2005.tif')

forest2010 <- raster('./Data/Forest - Vieilledent et al 2018/forest2010.tif')

forest2014 <- raster('./Data/Forest - Vieilledent et al 2018/forest2014.tif')

say('deforestation')

# convert NAs to 0

beginCluster(7)

zeroFunct <- function(x) ifelse(is.na(x), 0, x)

forest2000zeros <- clusterR(forest2000, calc, args=list(fun=zeroFunct))

forest2005zeros <- clusterR(forest2005, calc, args=list(fun=zeroFunct))

forest2010zeros <- clusterR(forest2010, calc, args=list(fun=zeroFunct))

forest2014zeros <- clusterR(forest2014, calc, args=list(fun=zeroFunct))

endCluster()

forest2000zeros <- humidForestBufferMask_utm38s * forest2000zeros

forest2005zeros <- humidForestBufferMask_utm38s * forest2005zeros

forest2010zeros <- humidForestBufferMask_utm38s * forest2010zeros

forest2014zeros <- humidForestBufferMask_utm38s * forest2014zeros

deforest2005 <- forest2000 - forest2005zeros

deforest2010 <- forest2005 - forest2010zeros

deforest2014 <- forest2010 - forest2014zeros

# convert 0s (deforested) to NAs

beginCluster(7)

fx <- function(x) ifelse(x == 1, 1, NA)

deforest2005 <- clusterR(deforest2005, calc, args=list(fun=fx))

deforest2010 <- clusterR(deforest2010, calc, args=list(fun=fx))

deforest2014 <- clusterR(deforest2014, calc, args=list(fun=fx))

endCluster()

names(deforest2005) <- 'deforest2005_utm38s'

names(deforest2010) <- 'deforest2010_utm38s'

names(deforest2014) <- 'deforest2014_utm38s'

writeRaster(deforest2005, './Data/Forest - Vieilledent et al 2018/deforest2005_utm38s', datatype='INT1U')

writeRaster(deforest2010, './Data/Forest - Vieilledent et al 2018/deforest2010_utm38s', datatype='INT1U')

writeRaster(deforest2014, './Data/Forest - Vieilledent et al 2018/deforest2014_utm38s', datatype='INT1U')

say('### distance to deforested patch ###')

say('2005')

distRast <- fasterRastDistance(deforest2005, metric='euclidean', meters=TRUE, fillNAs=TRUE, grassDir=grassDir)

distRast <- distRast * humidForestBufferMask_utm38s

distRast <- round(distRast)

names(distRast) <- 'distToDefoMeters2005_utm38s'

writeRaster(distRast, './Data/Forest - Vieilledent et al 2018/distToDefoMeters2005_utm38s', datatype='INT4S')

say('2010')

distRast <- fasterRastDistance(deforest2010, metric='euclidean', meters=TRUE, fillNAs=TRUE, grassDir=grassDir)

distRast <- distRast * humidForestBufferMask_utm38s

distRast <- round(distRast)

names(distRast) <- 'distToDefoMeters2010_utm38s'

writeRaster(distRast, './Data/Forest - Vieilledent et al 2018/distToDefoMeters2010_utm38s', datatype='INT4S')

say('2014')

distRast <- fasterRastDistance(deforest2014, metric='euclidean', meters=TRUE, fillNAs=TRUE, grassDir=grassDir)

distRast <- distRast * humidForestBufferMask_utm38s

distRast <- round(distRast)

names(distRast) <- 'distToDefoMeters2014_utm38s'

writeRaster(distRast, './Data/Forest - Vieilledent et al 2018/distToDefoMeters2014_utm38s', datatype='INT4S')

say('### distance to forest edge ###')

say('2000')

distRast <- fasterRastDistance(forest2000, metric='euclidean', meters=TRUE, fillNAs = FALSE, grassDir=grassDir)

distRast <- round(distRast)

distRast <- distRast * humidForestBufferMask_utm38s

names(distRast) <- 'distToForestEdgeMeters2000_utm38s'

writeRaster(distRast, './Data/Forest - Vieilledent et al 2018/distToForestEdgeMeters2000_utm38s', datatype='INT4S')

say('2005')

distRast <- fasterRastDistance(forest2005, metric='euclidean', meters=TRUE, fillNAs=FALSE, grassDir=grassDir)

distRast <- round(distRast)

distRast <- distRast * humidForestBufferMask_utm38s

names(distRast) <- 'distToForestEdgeMeters2005_utm38s'

writeRaster(distRast, './Data/Forest - Vieilledent et al 2018/distToForestEdgeMeters2005_utm38s', datatype='INT4S')

say('2010')

distRast <- fasterRastDistance(forest2010, metric='euclidean', meters=TRUE, fillNAs=FALSE, grassDir=grassDir)

distRast <- round(distRast)

distRast <- distRast * humidForestBufferMask_utm38s

names(distRast) <- 'distToForestEdgeMeters2010_utm38s'

writeRaster(distRast, './Data/Forest - Vieilledent et al 2018/distToForestEdgeMeters2010_utm38s', datatype='INT4S')

say('2014')

distRast <- fasterRastDistance(forest2014, metric='euclidean', meters=TRUE, fillNAs=FALSE, grassDir=grassDir)

distRast <- round(distRast)

distRast <- distRast * humidForestBufferMask_utm38s

names(distRast) <- 'distToForestEdgeMeters2014_utm38s'

writeRaster(distRast, './Data/Forest - Vieilledent et al 2018/distToForestEdgeMeters2014_utm38s', datatype='INT4S')

say('### forest fragmentation ###')

cores <- 4

say('2000')

frag2000 <- fasterFragmentation(rast = forest2000zeros, size = 5, pad = TRUE, padValue = NA, calcDensity = TRUE, calcConnect = TRUE, calcClass = TRUE, na.rm = TRUE, undet = 'perforated', cores = cores, forceMulti = TRUE)

frag2000 <- frag2000 * humidForestBufferMask_utm38s

frag2000[['density']] <- round(100 * frag2000[['density']])

frag2000[['connect']] <- round(100 * frag2000[['connect']])

names(frag2000) <- c('forestFragClass2000_utm38s', 'forestFragDensity2000_utm38s', 'forestFragConnect2000_utm38s')

writeRaster(frag2000[['forestFragClass2000_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragClass2000_utm38s', datatype='INT1U')

writeRaster(frag2000[['forestFragDensity2000_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragDensity2000_utm38s', datatype='INT1U')

writeRaster(frag2000[['forestFragConnect2000_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragConnect2000_utm38s', datatype='INT1U')

say('2005')

frag2005 <- fasterFragmentation(rast = forest2005zeros, size = 5, pad = TRUE, padValue = NA, calcDensity = TRUE, calcConnect = TRUE, calcClass = TRUE, na.rm = TRUE, undet = 'perforated', cores = cores, forceMulti = TRUE)

frag2005[['density']] <- round(100 * frag2005[['density']])

frag2005[['connect']] <- round(100 * frag2005[['connect']])

frag2005 <- frag2005 * humidForestBufferMask_utm38s

names(frag2005) <- c('forestFragClass2005_utm38s', 'forestFragDensity2005_utm38s', 'forestFragConnect2005_utm38s')

writeRaster(frag2005[['forestFragClass2005_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragClass2005_utm38s', datatype='INT1U')

writeRaster(frag2005[['forestFragDensity2005_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragDensity2005_utm38s', datatype='INT1U')

writeRaster(frag2005[['forestFragConnect2005_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragConnect2005_utm38s', datatype='INT1U')

say('2010')

frag2010 <- fasterFragmentation(rast = forest2010zeros, size = 5, pad = TRUE, padValue = NA, calcDensity = TRUE, calcConnect = TRUE, calcClass = TRUE, na.rm = TRUE, undet = 'perforated', cores = cores, forceMulti = TRUE)

frag2010[['density']] <- round(100 * frag2010[['density']])

frag2010[['connect']] <- round(100 * frag2010[['connect']])

frag2010 <- frag2010 * humidForestBufferMask_utm38s

names(frag2010) <- c('forestFragClass2010_utm38s', 'forestFragDensity2010_utm38s', 'forestFragConnect2010_utm38s')

writeRaster(frag2010[['forestFragClass2010_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragClass2010_utm38s', datatype='INT1U')

writeRaster(frag2010[['forestFragDensity2010_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragDensity2010_utm38s', datatype='INT1U')

writeRaster(frag2010[['forestFragConnect2010_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragConnect2010_utm38s', datatype='INT1U')

say('2014')

frag2014 <- fasterFragmentation(rast = forest2014zeros, size = 5, pad = TRUE, padValue = NA, calcDensity = TRUE, calcConnect = TRUE, calcClass = TRUE, na.rm = TRUE, undet = 'perforated', cores = cores, forceMulti = TRUE)

frag2014[['density']] <- round(100 * frag2014[['density']])

frag2014[['connect']] <- round(100 * frag2014[['connect']])

frag2014 <- frag2014 * humidForestBufferMask_utm38s

names(frag2014) <- c('forestFragClass2014_utm38s', 'forestFragDensity2014_utm38s', 'forestFragConnect2014_utm38s')

writeRaster(frag2014[['forestFragClass2014_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragClass2014_utm38s', datatype='INT1U')

writeRaster(frag2014[['forestFragDensity2014_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragDensity2014_utm38s', datatype='INT1U')

writeRaster(frag2014[['forestFragConnect2014_utm38s']], './Data/Forest - Vieilledent et al 2018/forestFragConnect2014_utm38s', datatype='INT1U')

say('### protected areas ###')

pa <- shapefile('./Data/Protected Areas/WDPA_Sept2018_MDG-shapefile-polygons')

# rasterize

pa_utm38s <- sp::spTransform(pa, CRS(madEaProj))

# crop to humid forest

humidForestBufferMask_utm38s <- raster('./Study Region & Masks/UTM 38S 30-m Resolution/humidForestBufferMask_utm38s.tif')

pa_utm38s <- crop(pa_utm38s, humidForestBufferMask_utm38s)

paRast_utm38s <- fasterRasterize(pa_utm38s, humidForestBufferMask_utm38s, grassDir=grassDir)

paRast_utm38s <- humidForestBufferMask_utm38s * paRast_utm38s

names(paRast_utm38s) <- 'protectedAreas_utm38s'

writeRaster(paRast_utm38s, './Data/Protected Areas/wdpa_utm38s')

# say('##################################################')

# say('### define areas of deforestation for modeling ###')

# say('##################################################')

# say('Divide forest/deforestation into geographic clusters using a clustering algorithm. Model each separately.')

# humidForestBufferMask_utm38s <- raster('./Study Region & Masks/UTM 38S 30-m Resolution/humidForestBufferMask_utm38s.tif')

# # generate sample of points from areas deforested 2000-2014

# forest2000 <- raster('./Data/Forest - Vieilledent et al 2018/forest2000.tif')

# forestSamples <- sampleRast(forest2000, n=10000, adjArea=FALSE, replace=FALSE, prob=FALSE)

# # cluster deforested sites

# forestSamples <- SpatialPoints(forestSamples, CRS(madEaProj))

# forestSamples <- sp::spTransform(forestSamples, getCRS('wgs84', TRUE))

# dists <- pointDist(forestSamples)

# dists <- as.dist(dists)

# library(cluster)

# groups <- agnes(dists)

# png('./Deforestation Models/Dendrogram of Forest Clusters Calculated Using UPGMA.png', width=900, height=900)

# plot(groups, main='Forest Clusters')

# dev.off()

# numFolds <- 6

# gFolds <- geoFold(forestSamples, numFolds)

# # assign cells to each group

# dirCreate('./Deforestation Models')

# for (i in 1:numFolds) {

# say('Calculating distance to fold ', i, '...')

# foldSamples <- defoSamples[gFolds == i]

# center <- gCentroid(foldSamples)

# distToFold <- fasterVectToRastDistance(humidForestBufferMask_utm38s, center)

# distToFold <- round(distToFold)

# names(distToFold) <- paste0('distToFold_', i)

# writeRaster(distToFold, paste0('./Deforestation Models/distToFold_', i), datatype='INT4S')

# }

# distToFolds <- raster::stack(listFiles('./Deforestation Models', pattern='distToFold'))

# beginCluster(7)

# foldAssign <- calc(distToFolds, which.min)

# endCluster()

# names(foldAssign) <- 'deforestationRegion'

# writeRaster(foldAssign, './Deforestation Models/deforestationRegion', datatype='INT1U')

# say('############################################################')

# say('### generate samples for modeling deforestation LOCATION ###')

# say('############################################################')

# say('Predictors:')

# say(' elevation')

# say(' slope')

# say(' forest fragmentation index')

# say(' distance to forest edge')

# say(' distance to previously deforested patch')

# say(' distance to road')

# say(' distance to settlement')

# say(' distance to nearest inland water body')

# say(' distance to nearest coast')

# say(' un/protected')

# say(' longitude and latitude')

# say(' human population density')

# # number of forested sites in starting time period... getting a little extra to account for NAs

# numTrainSites <- 40000

# numTestSites <- 10000

# numSites <- numTrainSites + numTestSites

# periodStarts <- c(2000, 2005, 2010)

# periodEnds <- c(2005, 2010, 2014)

# dirCreate('./Deforestation Models')

# ### by TIME PERIOD

# for (period in 1:3) {

# periodStart <- periodStarts[period]

# periodEnd <- periodEnds[period]

# say('Extracting deforestation predictors for ', periodStart, ' to ', periodEnd)

# ### generate sample sites and response variable (hard to place sites on the sparse forest rasters so getting more than needed)

# forestStart <- raster(paste0('./Data/Forest - Vieilledent et al 2018/forest', periodStart, '.tif'))

# forestEnd <- raster(paste0('./Data/Forest - Vieilledent et al 2018/forest', periodEnd, '.tif'))

# sites <- sampleRast(forestStart, 1.05 * numSites, adjArea=FALSE, replace=FALSE, prob=FALSE)

# # create response variable (0 = not deforested, 1 = deforested)

# coverStart <- extract(forestStart, sites)

# coverEnd <- extract(forestEnd, sites)

# defo <- rep(0, numTrainSites + numTestSites)

# defo[is.na(coverEnd)] <- 1

# ### extract predictors

# preds <- as.data.frame(sites)

# names(preds) <- c('longitude', 'latitude')

# # elevation

# rast <- raster('./Data/Topography - GMTED2010/elevationGmted2010_utm38s.tif')

# preds <- cbind(preds, extract(rast, sites))

# names(preds)[ncol(preds)] <- 'elevation'

# # slope

# rast <- raster('./Data/Topography - GMTED2010/slopeGmted2010_utm38s.tif')

# preds <- cbind(preds, extract(rast, sites))

# names(preds)[ncol(preds)] <- 'slope'

# # distance to nearest settlement

# rast <- raster('./Data/Geographic Places (NIMA from DIVA-GIS)/distToSettlementMeters_utm38s.tif')

# preds <- cbind(preds, extract(rast, sites))

# names(preds)[ncol(preds)] <- 'distToSettlement'

# # distance to nearest road

# rast <- raster('./Data/Roads/distToNearestRoadMeters_utm38s.tif')

# preds <- cbind(preds, extract(rast, sites))

# names(preds)[ncol(preds)] <- 'distToRoad'

# # distance to nearest inland water body

# rast <- raster('./Data/Water Bodies/distToNearestInlandWaterMeters_utm38s.tif')

# preds <- cbind(preds, extract(rast, sites))

# names(preds)[ncol(preds)] <- 'distToInlandWater'

# # distance to nearest coast

# rast <- raster('./Data/Water Bodies/distToNearestCoastMeters_utm38s.tif')

# preds <- cbind(preds, extract(rast, sites))

# names(preds)[ncol(preds)] <- 'distToCoast'

# # un/protected

# rast <- raster('./Data/Protected Areas/wdpa_utm38s.tif')

# preds <- cbind(preds, extract(rast, sites))

# names(preds)[ncol(preds)] <- 'protectedArea'

# preds$protectedArea[is.na(preds$protectedArea)] <- 0

# # forest fragmentation class

# rast <- raster(paste0('./Data/Forest - Vieilledent et al 2018/forestFragClass', periodStart, '_utm38s.tif'))

# preds <- cbind(preds, extract(rast, sites))

# names(preds)[ncol(preds)] <- 'fragClass'

# # distance to nearest deforested patch (not included for 2000-2005)

# if (periodStart != 2000) {

# rast <- raster(paste0('./Data/Forest - Vieilledent et al 2018/distToDefoMeters', periodStart, '_utm38s.tif'))

# preds <- cbind(preds, extract(rast, sites))

# names(preds)[ncol(preds)] <- 'distToDefo'

# } else {

# preds$distToDefo <- NA

# }

# # distance to nearest forest edge

# rast <- raster(paste0('./Data/Forest - Vieilledent et al 2018/distToForestEdgeMeters', periodStart, '_utm38s.tif'))

# preds <- cbind(preds, extract(rast, sites))

# names(preds)[ncol(preds)] <- 'distToForestEdge'

# # human population

# sites_utm38s <- SpatialPoints(sites, CRS(madEaProj))

# sites_wgs84 <- sp::spTransform(sites_utm38s, getCRS('wgs84', TRUE))

# rast <- raster(paste0('./Data/Population/popDensity', periodStart, '_wgs84.tif'))

# preds <- cbind(preds, extract(rast, sites_wgs84))

# names(preds)[ncol(preds)] <- 'popDensity'

# # collate data

# intervals <- data.frame(interval=rep(periodEnd - periodStart, nrow(preds)))

# thisData <- cbind(defo, intervals, preds)

# thisData$slope[is.na(thisData$slope)] <- 0

# nas <- naRows(thisData)

# if (length(nas) > 0) thisData <- thisData[1:numSites, ]

# thisData$startEndYear <- paste(periodStart, periodEnd)

# thisData$period <- period

# trainTest <- data.frame(trainTest = c(rep('train', numTrainSites), rep('test', numTestSites)))

# thisData <- as.data.frame(thisData)

# thisData <- cbind(trainTest, thisData)

# collatedData <- if (exists('collatedData')) {

# rbind(collatedData, thisData)

# } else {

# thisData

# }

# } # next time period

# collatedData$fragClass[collatedData$fragClass == 5] <- 4

# trainData <- collatedData[collatedData$trainTest == 'train', ]

# testData <- collatedData[collatedData$trainTest == 'test', ]

# save(trainData, file='./Deforestation Models/Training Sites & Predictors Located in Humid Forest UTM 38S.RData')

# save(testData, file='./Deforestation Models/Test Sites & Predictors Located in Humid Forest UTM 38S.RData')

# say('###############################################')

# say('### estimate deforestation VOLUME (vs rate) ###')

# say('###############################################')

# say('We are estimating deforestation VOLUME (vs rate) as the total area deforested each year as a function of total population size. Using rate assumes that the same *proportion* of forest is lost each year, whereas a static population should require the same amount of fiber product per year, regardless of stock. We will assume that the total amount of deforestation occurring in a year is a linear function of total population size such that: D = alpha * N0 + alpha * N1 + alpha * N2 + alpha * N3 + alpha * N4 = alpha * SUM N, where D is total area deforested between year Y0 and Y4, alpha is a constant, and N0, N1, N3, ... is the total population in years 0, 1, 2, ... for a given period.', breaks=80)

# load('./Data/Population Forecasts (UN)/WPP2017_PopulationBySingleAgeSex.RData')

# madPop <- worldPop[worldPop$Location == 'Madagascar', ]

# madPop <- summaryBy(PopTotal ~ Time, data=madPop, FUN=sum)

# names(madPop) <- c('year', 'population')

# deforest2005 <- raster('./Data/Forest - Vieilledent et al 2018/deforest2005_utm38s.tif')

# deforest2010 <- raster('./Data/Forest - Vieilledent et al 2018/deforest2010_utm38s.tif')

# deforest2014 <- raster('./Data/Forest - Vieilledent et al 2018/deforest2014_utm38s.tif')

# D1 <- cellStats(deforest2005, 'sum')

# D2 <- cellStats(deforest2010, 'sum')

# D3 <- cellStats(deforest2014, 'sum')

# # re-weight first and last years assuming population data was taken in middle of each year

# N1 <- madPop$population[madPop$year >= 2000 & madPop$year <= 2004]

# N1[1] <- 0.5 * N1[1]

# N1[length(N1)] <- 0.5 * N1[length(N1)]

# N2 <- madPop$population[madPop$year >= 2005 & madPop$year <= 2010]

# N2[1] <- 0.5 * N2[1]

# N2[length(N2)] <- 0.5 * N2[length(N2)]

# N3 <- madPop$population[madPop$year >= 2010 & madPop$year <= 2014]

# N3[1] <- 0.5 * N3[1]

# N3[length(N3)] <- 0.5 * N3[length(N3)]

# alpha1 <- D1 / sum(N1)

# alpha2 <- D2 / sum(N2)

# alpha3 <- D3 / sum(N3)

# ### predict forest amount loss per year... comes to more than total forest remaining in 2014, even if using *minimum* amount rate

# # MINIMUM

# defoAmount <- min(alpha1, alpha2, alpha3) * madPop$population[madPop$year >= 2015 & madPop$year <= 2080]

# defoAmount <- data.frame(year=2015:2080, amount_30x30m=defoAmount)

# dirCreate('./Figures & Tables/Forest - Deforestation Model Parameterization')

# write.csv(defoAmount, './Figures & Tables/Forest - Deforestation Model Parameterization/Deforested Amount per Year Using Minimum Loss Rate.csv', row.names=FALSE)

# ### predict forest amount loss per year... comes to more than total forest remaining in 2014, even if using *minimum* amount rate

# # MEAN

# defoAmount <- mean(alpha1, alpha2, alpha3) * madPop$population[madPop$year >= 2015 & madPop$year <= 2080]

# defoAmount <- data.frame(year=2015:2080, amount_30x30m=defoAmount)

# dirCreate('./Figures & Tables/Forest - Deforestation Model Parameterization')

# write.csv(defoAmount, './Figures & Tables/Forest - Deforestation Model Parameterization/Deforested Amount per Year Using Mean Loss Rate.csv', row.names=FALSE)

# ### predict forest amount loss per year... comes to more than total forest remaining in 2014, even if using *minimum* amount rate

# # MAXIMUM

# defoAmount <- max(alpha1, alpha2, alpha3) * madPop$population[madPop$year >= 2015 & madPop$year <= 2080]

# defoAmount <- data.frame(year=2015:2080, amount_30x30m=defoAmount)

# dirCreate('./Figures & Tables/Forest - Deforestation Model Parameterization')

# write.csv(defoAmount, './Figures & Tables/Forest - Deforestation Model Parameterization/Deforested Amount per Year Using Maximum Loss Rate.csv', row.names=FALSE)

# say('So... even if we use the minimum amount loss rate across the three periods, the total forest lost by 2080 is > all remaining forest.')

# say('##############################################')

# say('### Vielledent model of deforestation RATE ###')

# say('##############################################')

# load('./Deforestation Models/Training Sites & Predictors Located in Humid Forest UTM 38S.RData')

# model <- deforestation(defo ~ 1, interval=trainData$interval, data=trainData, burnin=8000, mcmc=2000, thin=1)

# mcmcMean <- as.matrix(model$mcmc)

# thetaMean <- mean(inv.logit(mcmcMean[ , 1])) # posterior mean of the annual deforestation rate

# dirCreate('./Figures & Tables/Forest - Deforestation Model Parameterization')

# write.csv(thetaMean, './Figures & Tables/Forest - Deforestation Model Parameterization/Mean Interval-Sensitive Deforestation Rate.csv', row.names=FALSE)

# say('##################################################')

# say('### Vielledent model of deforestation LOCATION ###')

# say('##################################################')

# load('./Deforestation Models/Training Sites & Predictors Located in Humid Forest UTM 38S.RData')

# # trainData <- trainData[!is.na(trainData$distToDefo), ]

# # transform

# trainData$longitude <- stretchMinMax(trainData$longitude, lower=-1, upper=1)

# trainData$latitude <- stretchMinMax(trainData$latitude, lower=-1, upper=1)

# trainData$distToSettlement <- log10(trainData$distToSettlement + 0.1)

# trainData$distToRoad <- log10(trainData$distToRoad + 0.1)

# trainData$distToCoast <- log10(trainData$distToCoast + 0.1)

# trainData$distToInlandWater <- log10(trainData$distToInlandWater + 0.1)

# trainData$distToDefo <- log10(trainData$distToDefo + 1)

# trainData$distToForestEdge <- log10(trainData$distToForestEdge + 0.1)

# variables <- c('longitude', 'latitude', 'elevation', 'slope', 'distToSettlement', 'distToRoad', 'distToInlandWater', 'distToCoast', 'as.factor(protectedArea)', 'as.factor(fragClass)', 'distToDefo', 'distToForestEdge', 'popDensity')

# trainData <- trainData[trainData$period != 1, ]

# ### evaluate predictor importance

# say('evaluate predictor importance')

# # formulas for all variables except one term

# forms <- list()

# for (v in seq_along(variables)) {

# forms[[v]] <- formula(paste('defo ~ ', paste(variables[-v], collapse='+'), sep=''))

# }

# # full model and null model

# forms[[length(forms) + 1]] <- formula(paste('defo ~ ', paste(variables, collapse='+'), sep=''))

# forms[[length(forms) + 1]] <- formula('defo ~ 1')

# # deviance

# dev <- rep(NA, length(forms))

# # train successive models

# for (i in seq_along(forms)) {

# say(forms[[i]])

# # train model

# model <- deforestation(forms[[i]], data=trainData, interval=trainData$interval, burnin=8000, mcmc=2000)

# mcmc <- as.matrix(model$mcmc)

# gammaHat <- apply(mcmc, 2 , mean)

# # predict to train data

# mfFixed <- model.frame(formula=forms[[i]], data=trainData)

# X <- model.matrix(attr(mfFixed, 'terms'), data=mfFixed)

# thetaHat <- inv.logit(X %*% gammaHat)

# thetaPrim <- 1 - (1 - thetaHat)^(trainData$interval)

# dev[i] <- -2 * sum(dbinom(trainData$defo, 1 , thetaPrim, 1))

# }

# # model comparison

# modelComp <- as.data.frame(matrix(NA, nrow=length(forms) + 1, ncol=2))

# names(modelComp) <- c('model', 'deviance')

# modelComp$model <- c(paste0('-', variables), 'full', 'null', 'saturated')

# modelComp$deviance <- c(round(dev, 3), 0)

# dirCreate('./Figures & Tables/Forest - Deforestation Model Parameterization')

# write.csv(modelComp, './Figures & Tables/Forest - Deforestation Model Parameterization/Deforestation Location Model Comparison.csv', row.names=FALSE)

# # variable importance (higher deviance ==> greater importance)

# varImp <- as.data.frame(matrix(NA, nrow=length(forms) - 1, ncol=3))

# names(varImp) <- c('model', 'deviance', 'varImp')

# varImp$model <- c(paste0('-', variables), 'full')

# varImp$deviance <- round(dev, 3)[1:(length(forms) - 1)]

# varImp$varImp <- (varImp$deviance - varImp$deviance[varImp$model == 'full'])

# write.csv(varImp, './Figures & Tables/Forest - Deforestation Model Parameterization/Deforestation Location Variable Importance.csv', row.names=FALSE)

# ### preliminary final model -- used to explore selected variables

# say('preliminary final model with selected variables')

# variables <- varImp$model[varImp$varImp > 0]

# variables <- substr(variables, 2, nchar(variables))

# form <- paste0('defo ~ ', paste0(variables, collapse=' + '))

# model <- deforestation(form, interval=trainData$interval, data=trainData, burnin=8000, mcmc=2000, verbose=1)

# summary(model$mcmc)

# plot(model$mcmc)

# ## full model with selected variables

# say('full model with selected variables')

# # removing distToCoast, distToInlandWater bc unexpected relationships to deforestation

# # removing fragmentation class bc correlated with distToForestEdge and because fragClass has non-zero coefficients only for frag class 4 (positive relationship)