[R-sig-Geo] Thresholds & reclassify raster

[John Wasige]

​Dear all​
,
I am requesting for your help to reclassify the gpp.tif (attached) raster
in R. Not sure my script here is doing the right thing.

​#Data structure​

class       : RasterBrick
dimensions  : 1243, 1409, 1751387, 1  (nrow, ncol, ncell, nlayers)
resolution  : 0.00225804, 0.00225804  (x, y)
extent      : -25.53282, -22.35125, 14.54232, 17.34906  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : L:\NDVI_WGS84\gpp.tif
names       :        gpp
min values  : -0.2176632
max values  :   2.989293


​###​

​I want to ​
use the statistical distribution of gpp.tif values set threshold for change
assignment. The thresholds should be set by the following two
​ ​
options:

(Option1) => use  mean and StD and assign all pixel values  < mean - StD as
negative change (class 1) and all values  > mean + StD as positive
change (class
3); all values in between would be assigned 'neutral' for the final
tabulation (class2).

(option2) => what may come to almost the same result is to do the
cumulative histogram of the calculated GPPproxy pixels and assign all
pixels below the 15 percentile to negative change (class1) and all pixel
above the 85 percentile to positive change (class3) and again the range in
between (15 to 85 percentile) would be assigned neutral (class2).

#R script

# classification
#Option1_mean_sd
gpp <- brick('L:/NDVI_WGS84/gpp.tif')
hist(as.matrix(GPP))
plot(GPPproxy)
m <- cellStats(gpp, 'mean')
sd <- cellStats(gpp, 'sd')
v <- m-sd
gppnew <- gpp-v
plot(gppnew)
gppnew[gppnew>0] <- 3
gppnew[gppnew<0] <- 1
gppnew[gppnew==0] <- 2
plot(gppnew)
writeRaster(gppnew, filename=paste("L:/NDVI_WGS84/gppnew_mean_sd",sep=''),
format="GTiff", overwrite=TRUE)

#Option2_percentile
max<- cellStats(gpp, 'max')
gppnew2 <-(gpp/max)*100
min <- cellStats(gpp, 'min')
m <- matrix(c(min, 15, 1,
              15, 85, 2,
              85, Inf, 3), ncol=3, byrow=TRUE)
gppnew_percentile <- reclassify(gppnew2, m)

hist(as.matrix(gppnew_percentile))
writeRaster(gppnew_percentile,
filename=paste("L:/NDVI_WGS84/gppnew_percentile",sep=''),
format="GTiff", overwrite=TRUE)
plot(gppnew_percentile)

​Thanks for your help

	[[alternative HTML version deleted]]