R Reprojecting global forest change rasters

Question

I have been trying to calculate forest cover and forest loss in West Kalimantan Province, Indonesia, during 2000-2020. I use the gfcanalysis package in R for this purpose. Unfortunately, the minimum value of the forest cover raster turned to be -20 (minus) after reprojecting it to UTM 49 S, whereas the original value should be 0 (zero).

Following http://azvoleff.com/articles/analyzing-forest-change-with-gfcanalysis/, here are my codes to generate the output:

library('sp')
library('raster')
library('rgdal')
library('gfcanalysis')
// set working directory
setwd("D:/Users/GFC")
getwd()
// prepare output folder
output_folder <- "C:/Users/GFC/forestloss"
// download boundary of West Papua Province (WPP) via GADM 
idn <- getData('GADM', country='IDN', level=1, download = TRUE)
// inspect GADM and choose WPP
idn@data
idn <- idn[idn$NAME_1 == "Kalimantan Barat",]
// reproject GADM to UTM 49S
idn.utm <- spTransform(idn, CRS("+init=epsg:32749"))
plot(idn.utm)
// calculate tiles needed to cover WPP
tiles <- calc_gfc_tiles(idn)
print(length(tiles))
// download GFC for WPP
download_tiles(tiles, output_folder)
// extract GFC data for WPP
gfc_extract <- extract_gfc(idn, output_folder, to_UTM=F, filename="GFC_extract.tif")
// reproject GFC
gfc_extract_utm <- projectRaster(gfc_extract, crs = "+init=epsg:32749")
gfc_extract_utm
// set forest threshold: 30%
forest_threshold <- 30
// thresholded GFC 
gfc_ extract_utm_th<- threshold_gfc(gfc_extract_utm, forest_threshold=forest_threshold, filename="gfc_ extract_utm_th.tif")
// calculate forest statistics 
gfc_statistics <- gfc_stats(idn.utm, gfc_ extract_utm_th)

Any idea what is happening?

Answer 1

Reprojection involves interpolation- this almost always results in your destination values being different from your source values. This is explained well here and here.

You need to take care to select an interpolation method that is appropriate for the type of data that you are working with. Since the default method for projectRaster() is method='bilinear' (appropriate for continuous variables), you might try method='ngb' (nearest neighbor; appropriate for categorical variables) instead and compare your result.

Below is a brief illustration:

#example data
r <- raster::raster(nrows=36, 
                    ncols=36, 
                    xmn=-121.6000, 
                    xmx=-121.5900, 
                    ymn=50.47000,
                    ymx=50.48000, 
                    resolution=0.0002778, 
                    crs = CRS('+proj=longlat +datum=WGS84 +no_defs'))
#set up some dummy values
r[]<-1:ncell(r)
r
class      : RasterLayer
dimensions : 36, 36, 1296  (nrow, ncol, ncell)
resolution : 0.0002778, 0.0002778  (x, y)
extent     : -121.6, -121.59, 50.47, 50.48  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
source     : memory
names      : layer
values     : 1, 1296  (min, max)
#bilinear (default) interpolation - appropriate for continuous variables
b<-projectRaster(r, method="bilinear",crs = CRS('+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0'))
b
class      : RasterLayer
dimensions : 45, 51, 2295  (nrow, ncol, ncell)
resolution : 19.7, 30.9  (x, y)
extent     : 173504.8, 174509.5, 5600849, 5602239  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs
source     : memory
names      : layer
values     : -4.586387, 1291.587  (min, max)
#nearest neighbor interpolation - appropriate for categorical variables
n<-projectRaster(r, method="ngb",crs = CRS('+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0'))
n
class      : RasterLayer
dimensions : 45, 51, 2295  (nrow, ncol, ncell)
resolution : 19.7, 30.9  (x, y)
extent     : 173504.8, 174509.5, 5600849, 5602239  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs
source     : memory
names      : layer
values     : 1, 1296  (min, max)

Use caution though- it may be more appropriate to have max and min values that are "off" rather than using the wrong interpolation method for your type of data!