[R-sig-Geo] Thresholds & reclassify raster

[Ben Tupper]

Hi,

Your attachment didn't make it through.

I'm not sure about your first option, but your second can use a combination of quantile() and findInterval() to create a classified raster.  Perhaps like this...


library(raster)
r <- raster()
r <- setValues(r, runif(ncell(r), 40, 90))

pp <- quantile(r, c(0, 0.15, 0.85))
# pp
#       0%      15%      85% 
# 40.00199 47.64569 82.50751 

ix <- findInterval(getValues(r), pp)

classified_r <- setValues(r, ix)

# > r
# class       : RasterLayer 
# dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
# resolution  : 1, 1  (x, y)
# extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory
# names       : layer 
# values      : 40.00199, 89.99892  (min, max)

# > classified_r
# class       : RasterLayer 
# dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
# resolution  : 1, 1  (x, y)
# extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory
# names       : layer 
# values      : 1, 3  (min, max)


Cheers,
Ben
> On Jun 22, 2017, at 4:04 AM, John Wasige <johnwasige at gmail.com> wrote:
> 
> ​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]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org