Thresholds & reclassify raster

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: 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 <[hidden email]> 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
> [hidden email]
> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Thresholds & reclassify raster

John Wasige
Many thanks Ben Tupper!
Best Rgds John

On Thu, Jun 22, 2017 at 1:47 PM, Ben Tupper <[hidden email]> 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 <[hidden email]> 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
> > [hidden email]
> > 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Loading...