[R-sig-Geo] Thresholds & reclassify raster

[John Wasige]

Many thanks Ben Tupper!
Best Rgds John

On Thu, Jun 22, 2017 at 1:47 PM, Ben Tupper <btupper at bigelow.org> wrote:

> 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
>
>
>
>


-- 
John Wasige
"There are no REGRATES in LIFE, just lessons (Jennifer Aniston)”

	[[alternative HTML version deleted]]