0

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?

2
  • Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
    – Community Bot
    Sep 14 at 16:03
  • Thanks for letting me know. I have revised my question. Hopefully, it's clear now.
    – Ing
    Sep 14 at 16:38

1 Answer 1

0

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!

1
  • Thanks for your suggestion, Roman. I obtain comparable results after using the nearest neighbor approach as the interpolation method.
    – Ing
    Sep 15 at 4:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.