Approach

Study site

alt text

alt text

Plotting the magnitude

## Loading required package: sp
## rgdal: version: 0.9-1, (SVN revision 518)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
## Path to GDAL shared files: /usr/share/gdal/1.10
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)
## 
## Attaching package: 'bfast'
## 
## The following object is masked _by_ '.GlobalEnv':
## 
##     dates
## 
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Loading required package: ncdf4
## Loading required package: maps
## Loading required package: mapdata
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/christopher/Documents/Master Thesis/PRODES/Manicore", layer: "clip_clearcut"
## with 15847 features and 13 fields
## Feature type: wkbPolygon with 2 dimensions
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/christopher/Documents/Master Thesis/PRODES/Manicore", layer: "clip_clearcut_greater_equals_2011_wgs"
## with 2337 features and 13 fields
## Feature type: wkbPolygon with 2 dimensions

Confustion matrix / error matrix

# analyze only the negative magnitudes because the EVI timeseries is a negative deveation from the model 
# all values become 1 between -0.6 and 0 and all values become 0 between 0 and 0.8
manicoreMagRasterRecl = reclassify(manicoreMagRaster, rcl=c(-0.6, 0, 1, 0, 0.8, 0))

# Binary map showing all negative magnitudes as green and all positve magnitueds as red
plot(manicoreMagRasterRecl, main="Green: all negative magnitudes \n White: all positive magnitudes")
plot(prodesPolygon, add=T, lwd=0.5)
plot(prodesPolygonSince2011, border="red", add=T, lwd=0.5)

#resetting the raster values
dummyRaster = reclassify(manicoreMagRasterRecl, rcl=c(-10000, 10000, 0))

# rasterize the prodes polygon where deforestation was monitored between 2011 and 2013
# update the value of the cells that overlap with the prodes polygon with the value 1
prodesPolygonSince2011Raster = rasterize(prodesPolygonSince2011, dummyRaster, update=TRUE, field=1)
plot(prodesPolygonSince2011Raster, main="Green: where deforestaion was \n monitored between 2011 and 2013")

prodes = getValues(prodesPolygonSince2011Raster) # prodes
bfm = getValues(manicoreMagRasterRecl) # bfm

confusionMatrix = aggregate(x=prodes, by=list(prodes, bfm), FUN = length)

confusionMatrix
##   Group.1 Group.2     x
## 1       0       0 29075
## 2       1       0  1053
## 3       0       1 30102
## 4       1       1  1330

Prodes is in the horizontal row and bfastmonitor is in the vertical column.

PRODES vs bfm 0 1 Row total
0 29068 1060 30128
1 30099 1333 31432
column total 59167 2393 61560

P(0|0) + P(1|1) = 30401

The overall accuracy for the negative magnitudes:

\(30401 / 61560 = 0.49\)

Time series

options(scipen=999)
prodesMaskSince2011 = mask(manicoreMagRaster, prodesPolygonSince2011)
plot(prodesMaskSince2011)

magVals = na.omit(getValues(prodesMaskSince2011))

plot(prodesMaskSince2011, col = rainbow(5))

fj5 = classIntervals(magVals, n=5, style="fisher", dataPrecision=2)
plot(prodesMaskSince2011, col = rainbow(5), breaks=round(fj5$brks, 2))
plot(prodesPolygonSince2011, add=TRUE)

# negative magnitudes
oranges <- brewer.pal(n=5, name="Oranges")

# http://colorbrewer2.org/
redPurples1 = c("#feebe2", "#fbb4b9", "#f768a1", "#c51b8a", "#7a0177")
redPurples2 = c("#f1eef6", "#d7b5d8", "#df65b0", "#dd1c77", "#980043")
redPurples3 = c("#c994c7", "#df65b0", "#e7298a", "#ce1256", "#980043")

fj5NegMag = classIntervals(magVals[magVals < 0], n=5, style="fisher", dataPrecision=2)
plot(prodesMaskSince2011, col = rev(redPurples3), zlim=c(round(fj5NegMag$brks[1], 2), 0), breaks=round(fj5NegMag$brks, 2),
     main="Negative magnitudes of bfastmonitor \n inside PRODES polygons from 2011 to 2013")
plot(prodesPolygonSince2011, add=TRUE)

# positive magnitudes
greens <- brewer.pal(n=5, name="Greens")
yellowGreens = c("#ffffcc", "#d9f0a3", "#addd8e", "#78c679", "#41ab5d")
yellowGreens2 = c("#addd8e", "#78c679", "#41ab5d", "#238443", "#005a32")

fj5PosMag = classIntervals(magVals[magVals > 0], n=5, style="fisher", dataPrecision=2)
plot(prodesMaskSince2011, col = yellowGreens2, zlim=c(0, round(fj5PosMag$brks[5], 2)), breaks=round(fj5PosMag$brks, 2),
     main="Positive magnitudes of bfastmonitor \n inside PRODES polygons from 2011 to 2013")
plot(prodesPolygonSince2011, add=TRUE)

# positve and negativ magnitudes
plot(prodesMaskSince2011, col = rev(redPurples3), zlim=c(round(fj5NegMag$brks[1], 2), 0), breaks=round(fj5NegMag$brks, 2),
    main="Positive (green) and negative (purple) magnitudes \n of bfastmonitor inside PRODES polygons from 2011 to 2013",
    legend=FALSE)
plot(prodesMaskSince2011, col = yellowGreens2, zlim=c(0, round(fj5PosMag$brks[5], 2)), breaks=round(fj5PosMag$brks, 2), legend=FALSE,
     add=TRUE)
plot(prodesPolygonSince2011, add=TRUE)

Time Series

############################
## Look at time series
############################

rasterstackcrop = "/home/christopher/Documents/git/MasterThesis/spatialAccuracy/manicore.gri"

rasterStack = stack(rasterstackcrop)

# scale the EVI values
rasterStack = calc(rasterStack, function(x) x/10000)

# get layer names from raster stack
layerNames = names(rasterStack)

# extracting julian date of layer name 
dates = lapply(layerNames, function(x) str_extract(x, perl("(?<=X).*?(?=.MOD|$)")))

# transform julian date to date
dates = as.Date(decipher.M3.date(as.numeric(unlist(dates))))

rasterStackTs = rts(rasterStack, dates)

####################
# Analyzing magnitudes inside PRODES polygons from 2011-2013
####################

####################
# analyzing the maximum magnitude value inside PRODES polygon from 2011-2013
####################

plot(prodesMaskSince2011, col = rev(redPurples3), zlim=c(round(fj5NegMag$brks[1], 2), 0), breaks=round(fj5NegMag$brks, 2),
    main="Location maximum magnitude value inside PRODES polygon",
    legend=FALSE)

plot(prodesMaskSince2011, col = yellowGreens2, zlim=c(0, round(fj5PosMag$brks[5], 2)), breaks=round(fj5PosMag$brks, 2), legend=FALSE,
     add=TRUE)

maxMagCell = which.max(prodesMaskSince2011)
maxMagCoords = xyFromCell(prodesMaskSince2011, maxMagCell)
getValues(prodesMaskSince2011)[which.max(prodesMaskSince2011)]
## [1] 0.5524009
points(maxMagCoords, col="yellow", lwd=1.5, cex=1.5, pch=1)

timeSeries = ts(as.numeric(rasterStackTs[maxMagCell,]), start=c(2000, 4), end=c(2014, 17), frequency=23)  
plot(timeSeries)

bfm = bfastmonitor(timeSeries, start=c(2010, 23), history = c("ROC"))
bfm
## 
## BFAST monitoring
## 
## 1. History period
## Stable period selected: 2009(11)--2010(22)
## Length (in years): 1.478261
## Model fit:
##   (Intercept)         trend    harmoncos1    harmoncos2    harmoncos3 
##  2.4285710512 -0.0083438309  0.0278285381  0.0199557342  0.0201784541 
##    harmonsin1    harmonsin2    harmonsin3 
##  0.0180200867  0.0003369443  0.0026576597 
## R-squared: 0.791965
## 
## 
## 2. Monitoring period
## Monitoring period assessed: 2010(23)--2014(17)
## Length (in years): 3.739130
## Break detected at: 2011(10)
plot(bfm)

####################
# analyzing the minumum magnitude value inside PRODES poylgons from 2011-2013
####################

plot(prodesMaskSince2011, col = rev(redPurples3), zlim=c(round(fj5NegMag$brks[1], 2), 0), breaks=round(fj5NegMag$brks, 2),
    main="Location minimum magnitude value inside PRODES polygon", legend=FALSE)

plot(prodesMaskSince2011, col = yellowGreens2, zlim=c(0, round(fj5PosMag$brks[5], 2)), breaks=round(fj5PosMag$brks, 2), legend=FALSE,
     add=TRUE)

minMagCell = which.min(prodesMaskSince2011)
minMagCoords = xyFromCell(prodesMaskSince2011, minMagCell)
getValues(prodesMaskSince2011)[which.min(prodesMaskSince2011)]
## [1] -0.2721663
points(minMagCoords, col="blue", lwd=1, cex=.5, pch=1)

timeSeries = ts(as.numeric(rasterStackTs[minMagCell,]), start=c(2000, 4), end=c(2014, 17), frequency=23)  
plot(timeSeries)

bfm = bfastmonitor(timeSeries, start=c(2010, 23), history = c("ROC"))
bfm
## 
## BFAST monitoring
## 
## 1. History period
## Stable period selected: 2009(3)--2010(22)
## Length (in years): 1.826087
## Model fit:
##  (Intercept)        trend   harmoncos1   harmoncos2   harmoncos3 
## -0.339224386  0.003825738 -0.021441499 -0.021532746 -0.009238024 
##   harmonsin1   harmonsin2   harmonsin3 
##  0.084511224  0.002629814  0.006166541 
## R-squared: 0.823723
## 
## 
## 2. Monitoring period
## Monitoring period assessed: 2010(23)--2014(17)
## Length (in years): 3.739130
## Break detected at: 2011(5)
plot(bfm)


#####################
# get the cell number to analyze
#####################

click(prodesMaskSince2011, cell=T)

####################
# pixel of large green patch inside a PRODES polygons from 2011-2013
####################

# the cell 44905 large green patch

plot(prodesMaskSince2011, col = rev(redPurples3), zlim=c(round(fj5NegMag$brks[1], 2), 0), breaks=round(fj5NegMag$brks, 2),
    main="Location minimum magnitude value inside PRODES polygon", legend=FALSE)

plot(prodesMaskSince2011, col = yellowGreens2, zlim=c(0, round(fj5PosMag$brks[5], 2)), breaks=round(fj5PosMag$brks, 2), legend=FALSE,
     add=TRUE)

cellOfGreenPatch = xyFromCell(prodesMaskSince2011, 44905)
points(cellOfGreenPatch, col="blue", lwd=1, cex=.5, pch=1)

getValues(prodesMaskSince2011)[44905]
## [1] 0.116781
timeSeries = ts(as.numeric(rasterStackTs[44905,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
plot(timeSeries)

bfm = bfastmonitor(timeSeries, start=c(2010, 23), history = c("ROC"))
bfm
## 
## BFAST monitoring
## 
## 1. History period
## Stable period selected: 2008(13)--2010(22)
## Length (in years): 2.391304
## Model fit:
##   (Intercept)         trend    harmoncos1    harmoncos2    harmoncos3 
##  0.8203756671 -0.0010681861  0.0724687961  0.0005235339  0.0022189069 
##    harmonsin1    harmonsin2    harmonsin3 
##  0.0581462951  0.0217831484  0.0072085843 
## R-squared: 0.891425
## 
## 
## 2. Monitoring period
## Monitoring period assessed: 2010(23)--2014(17)
## Length (in years): 3.739130
## Break detected at: 2011(2)
plot(bfm)

####################
# analyze magnitudes outside PRODES polygons (the inverse mask)
####################

prodesMaskSince2011Inv = mask(manicoreMagRaster, prodesPolygonSince2011, inverse=TRUE)
plot(prodesMaskSince2011Inv)

magValsInv = na.omit(getValues(prodesMaskSince2011Inv))

fj5 = classIntervals(magValsInv, n=5, style="fisher", dataPrecision=2)
plot(prodesMaskSince2011Inv, col = rainbow(5), breaks=round(fj5$brks, 2), main="Magnitudes outside of PRODES polygons 2011-2013")
plot(prodesPolygonSince2011, add=TRUE)

# http://colorbrewer2.org/
redPurples1 = c("#feebe2", "#fbb4b9", "#f768a1", "#c51b8a", "#7a0177")
redPurples2 = c("#f1eef6", "#d7b5d8", "#df65b0", "#dd1c77", "#980043")
redPurples3 = c("#c994c7", "#df65b0", "#e7298a", "#ce1256", "#980043")

fj5NegMagInv = classIntervals(magValsInv[magValsInv < 0], n=5, style="fisher", dataPrecision=2)
plot(prodesMaskSince2011Inv, col = rev(redPurples3), zlim=c(round(fj5NegMagInv$brks[1], 2), 0), breaks=round(fj5NegMagInv$brks, 2),
     main="Negative magnitudes of bfastmonitor \n outside PRODES polygons from 2011 to 2013")
plot(prodesPolygonSince2011, add=TRUE)

# positive magnitudes
greens <- brewer.pal(n=5, name="Greens")
yellowGreens = c("#ffffcc", "#d9f0a3", "#addd8e", "#78c679", "#41ab5d")
yellowGreens2 = c("#addd8e", "#78c679", "#41ab5d", "#238443", "#005a32")

fj5PosMagInv = classIntervals(magValsInv[magValsInv > 0], n=5, style="fisher", dataPrecision=2)
plot(prodesMaskSince2011Inv, col = yellowGreens2, zlim=c(0, round(fj5PosMagInv$brks[5], 2)), breaks=round(fj5PosMagInv$brks, 2),
     main="Positive magnitudes of bfastmonitor \n outside PRODES polygons from 2011 to 2013")
plot(prodesPolygonSince2011, add=TRUE)

# positve and negativ magnitudes
plot(prodesMaskSince2011Inv, col = rev(redPurples3), zlim=c(round(fj5NegMagInv$brks[1], 2), 0), breaks=round(fj5NegMagInv$brks, 2),
     main="Positive and negative magnitudes of bfastmonitor \n outside PRODES polygons from 2011 to 2013", legend=FALSE)
plot(prodesMaskSince2011Inv, col = yellowGreens2, zlim=c(0, round(fj5PosMagInv$brks[5], 2)), breaks=round(fj5PosMagInv$brks, 2),
     legend=FALSE, add=TRUE)
plot(prodesPolygon, border="blue", lwd=0.5, add=TRUE)
plot(prodesPolygonSince2011, add=TRUE)



#################
# pixel of large green patch outside PRODES polygons 2011-2013 but inside PRODES polygon reported before 2011
#################

click(prodesMaskSince2011Inv, cell=TRUE)

plot(prodesMaskSince2011Inv, col = rev(redPurples3), zlim=c(round(fj5NegMagInv$brks[1], 2), 0), breaks=round(fj5NegMagInv$brks, 2),
     main="Positive and negative magnitudes of bfastmonitor \n outside PRODES polygons from 2011 to 2013", legend=FALSE)
plot(prodesMaskSince2011Inv, col = yellowGreens2, zlim=c(0, round(fj5PosMagInv$brks[5], 2)), breaks=round(fj5PosMagInv$brks, 2),
     legend=FALSE, add=TRUE)

cellOfGreenPatch = xyFromCell(prodesMaskSince2011Inv, 24125)
points(cellOfGreenPatch, col="blue", lwd=1, cex=.5, pch=1)

getValues(prodesMaskSince2011Inv)[24125]
## [1] 0.153106
timeSeries = ts(as.numeric(rasterStackTs[24125,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
plot(timeSeries)

bfm = bfastmonitor(timeSeries, start=c(2010, 23), history = c("ROC"))
bfm
## 
## BFAST monitoring
## 
## 1. History period
## Stable period selected: 2008(17)--2010(22)
## Length (in years): 2.217391
## Model fit:
##   (Intercept)         trend    harmoncos1    harmoncos2    harmoncos3 
##  1.0502913201 -0.0024785620  0.0593529984  0.0252885978  0.0037108760 
##    harmonsin1    harmonsin2    harmonsin3 
##  0.1089425767  0.0158039834  0.0002679368 
## R-squared: 0.838737
## 
## 
## 2. Monitoring period
## Monitoring period assessed: 2010(23)--2014(17)
## Length (in years): 3.739130
## Break detected at: 2011(3)
plot(bfm)

#################
# pixel of large purple patch outside PRODES polygons 2011-2013 but inside PRODES polygon reported before 2011
#################

plot(prodesMaskSince2011Inv, col = rev(redPurples3), zlim=c(round(fj5NegMagInv$brks[1], 2), 0), breaks=round(fj5NegMagInv$brks, 2),
     main="Positive and negative magnitudes of bfastmonitor \n outside PRODES polygons from 2011 to 2013", legend=FALSE)
plot(prodesMaskSince2011Inv, col = yellowGreens2, zlim=c(0, round(fj5PosMagInv$brks[5], 2)), breaks=round(fj5PosMagInv$brks, 2),
     legend=FALSE, add=TRUE)

cellOfPurplePatch = xyFromCell(prodesMaskSince2011Inv, 43644 )
points(cellOfPurplePatch, col="blue", lwd=1, cex=.5, pch=1)

getValues(prodesMaskSince2011Inv)[43644 ]
## [1] -0.147182
timeSeries = ts(as.numeric(rasterStackTs[43644,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
plot(timeSeries)

bfm = bfastmonitor(timeSeries, start=c(2010, 23), history = c("ROC"))
bfm
## 
## BFAST monitoring
## 
## 1. History period
## Stable period selected: 2006(21)--2010(22)
## Length (in years): 4.043478
## Model fit:
##   (Intercept)         trend    harmoncos1    harmoncos2    harmoncos3 
##  0.5048782948  0.0002735842  0.0224848973 -0.0063351709 -0.0036853405 
##    harmonsin1    harmonsin2    harmonsin3 
##  0.0657704571  0.0051907348  0.0052146720 
## R-squared: 0.318991
## 
## 
## 2. Monitoring period
## Monitoring period assessed: 2010(23)--2014(17)
## Length (in years): 3.739130
## Break detected at: 2011(21)
plot(bfm)

#################
# pixel of large green patch outside of any PRODES polygons 2000-2013
#################

plot(prodesMaskSince2011Inv, col = rev(redPurples3), zlim=c(round(fj5NegMagInv$brks[1], 2), 0), breaks=round(fj5NegMagInv$brks, 2),
     main="Positive and negative magnitudes of bfastmonitor \n outside PRODES polygons from 2011 to 2013", legend=FALSE)
plot(prodesMaskSince2011Inv, col = yellowGreens2, zlim=c(0, round(fj5PosMagInv$brks[5], 2)), breaks=round(fj5PosMagInv$brks, 2),
     legend=FALSE, add=TRUE)
plot(prodesPolygon, border="blue", lwd=0.5, add=TRUE)
plot(prodesPolygonSince2011, add=TRUE)

Magnitude Type inside PRODES 2000-2013 outside PRODES 2000-2013 inside PRODES 2011-2013 outside PRODES 2011-2013
positive Mag Content Cell Content Cell Time series indicates shows that deforestation happenend Content Cell
negative Mag Content Cell Content Cell Content Cell Content Cell