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 <johnwas...@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@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 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo