Count values "greater than" in rasterbrick layers

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

Count values "greater than" in rasterbrick layers

Thiago V. dos Santos
Dear all,

I have a 31-layer rasterbrick where each layer represents daily rainfall records for a certain month.

My goal is to perform two calculations in this brick:

a) Sum up total rainfall for the entire month

b) For the entire month (all rasterbrick layers), calculate the number of days where it rained (i.e. rainfall > 1)

Problem a) is straightforward and I solved it with the help of a simple calc command. 

Problem b), however, is where I am having troubles. The command I came up with doesn't work, as the final raster has values way greater than 31 (which would be the limit of each cell, if it had rained every day of the month on that cell).

This is my code so far. Any suggestions?

library(raster)

# generate raster data
r <- raster(ncol=384, nrow=190)
b <- brick(sapply(1:31, function(i) setValues(r, rnorm(ncell(r), i, 3))))

# sum up rain
acc.rain <- calc(b, fun=sum)

# count wet days
wet.days <- calc(b, fun=function(x,na.rm) sum(x[x>1]))
 

Many thanks in advance,
--
Thiago V. dos Santos
PhD student
Land and Atmospheric Science
University of Minnesota
http://www.laas.umn.edu/CurrentStudents/MeettheStudents/ThiagodosSantos/index.htm
Phone: (612) 323 9898
        [[alternative HTML version deleted]]


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

Re: Count values "greater than" in rasterbrick layers

Stachelek, Joseph
Hi Thiago,

Ignoring the fact that you have negative rainfall amounts in your example, I would go about this problem by reclassifying the raster brick to a binary classification. Then if you sum you get raw counts.

#count wet days
b.wet<-reclassify(b,c(-Inf,0,0))
b.wet<-reclassify(b.wet,c(0,+Inf,1))
plot(b.wet[[1]])
plot(sum(b.wet))

-----Original Message-----
From: [hidden email] [mailto:[hidden email]] On Behalf Of Thiago V. dos Santos
Sent: Wednesday, July 16, 2014 5:51 AM
To: R-SIG list
Subject: [R-sig-Geo] Count values "greater than" in rasterbrick layers

Dear all,

I have a 31-layer rasterbrick where each layer represents daily rainfall records for a certain month.

My goal is to perform two calculations in this brick:

a) Sum up total rainfall for the entire month

b) For the entire month (all rasterbrick layers), calculate the number of days where it rained (i.e. rainfall > 1)

Problem a) is straightforward and I solved it with the help of a simple calc command.

Problem b), however, is where I am having troubles. The command I came up with doesn't work, as the final raster has values way greater than 31 (which would be the limit of each cell, if it had rained every day of the month on that cell).

This is my code so far. Any suggestions?

library(raster)

# generate raster data
r <- raster(ncol=384, nrow=190)
b <- brick(sapply(1:31, function(i) setValues(r, rnorm(ncell(r), i, 3))))

# sum up rain
acc.rain <- calc(b, fun=sum)

# count wet days
wet.days <- calc(b, fun=function(x,na.rm) sum(x[x>1]))


Many thanks in advance,
--
Thiago V. dos Santos
PhD student
Land and Atmospheric Science
University of Minnesota
http://www.laas.umn.edu/CurrentStudents/MeettheStudents/ThiagodosSantos/index.htm
Phone: (612) 323 9898
        [[alternative HTML version deleted]]



We value your opinion. Please take a few minutes to share your comments on the service you received from the District by clicking on this link<http://my.sfwmd.gov/portal/page/portal/pg_grp_surveysystem/survey%20ext?pid=1653>.

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

Re: Count values "greater than" in rasterbrick layers

Thiago V. dos Santos
Hi Joseph,

Thanks a lot for your tip - it completely solved my problem.
 
Greetings,
--
Thiago V. dos Santos
PhD student
Land and Atmospheric Science
University of Minnesota
http://www.laas.umn.edu/CurrentStudents/MeettheStudents/ThiagodosSantos/index.htm
Phone: (612) 323 9898


On Wednesday, July 16, 2014 7:30 AM, "Stachelek, Joseph" <[hidden email]> wrote:
 


Hi Thiago,

Ignoring the fact that you have negative rainfall amounts in your example, I would go about this problem by reclassifying the raster brick to a binary classification. Then if you sum you get raw counts.

#count wet days
b.wet<-reclassify(b,c(-Inf,0,0))
b.wet<-reclassify(b.wet,c(0,+Inf,1))
plot(b.wet[[1]])
plot(sum(b.wet))


-----Original Message-----
From: [hidden email] [mailto:[hidden email]] On Behalf Of Thiago V. dos Santos
Sent: Wednesday, July 16, 2014 5:51 AM
To: R-SIG list
Subject: [R-sig-Geo] Count values "greater than" in rasterbrick layers

Dear all,

I have a 31-layer rasterbrick where each layer represents daily rainfall records for a certain month.

My goal is to perform two calculations in this brick:

a) Sum up total rainfall for the entire month

b) For the entire month (all rasterbrick layers), calculate the number of days where it rained (i.e. rainfall > 1)

Problem a) is straightforward and I solved it with the help of a simple calc command.

Problem b), however, is where I am having troubles. The command I came up with doesn't work, as the final raster has values way greater than 31 (which would be the limit of each cell, if it had rained every day of the month on that cell).

This is my code so far. Any suggestions?

library(raster)

# generate raster data
r <- raster(ncol=384, nrow=190)
b <- brick(sapply(1:31, function(i) setValues(r, rnorm(ncell(r), i, 3))))

# sum up rain
acc.rain <- calc(b, fun=sum)

# count wet days
wet.days <- calc(b, fun=function(x,na.rm) sum(x[x>1]))


Many thanks in advance,
--
Thiago V. dos Santos
PhD student
Land and Atmospheric Science
University of Minnesota
http://www.laas.umn.edu/CurrentStudents/MeettheStudents/ThiagodosSantos/index.htm
Phone: (612) 323 9898
        [[alternative HTML version deleted]]



We value your opinion. Please take a few minutes to share your comments on the service you received from the District by clicking on this link<http://my.sfwmd.gov/portal/page/portal/pg_grp_surveysystem/survey%20ext?pid=1653
>.
        [[alternative HTML version deleted]]


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

Re: Count values "greater than" in rasterbrick layers

dschneiderch
In reply to this post by Thiago V. dos Santos
Thiago -
With regards to your original attempt for part b, I  believe using length() instead of sum() will get you what you want.  In this case, sum(x) and sum(x[x>1]) are the same assuming you have no values less than 1. Length can be used as a proxy for counting how many days rain>1.
Reply | Threaded
Open this post in threaded view
|

Re: Count values "greater than" in rasterbrick layers

Robert Hijmans
Thiago,
Here are more direct approaches:

# generate raster data
r <- raster(ncol=384, nrow=190)
b <- brick(sapply(1:31, function(i) setValues(r, rnorm(ncell(r), i, 3))))

# sum up rain
acc.rain <- sum(b)

# count wet days
wet.days <- sum(b > 0)

or with calc
wet.days <-  calc(b, function(x) sum(x > 0))

Robert

On Tue, Jul 22, 2014 at 1:36 PM, dschneiderch
<[hidden email]> wrote:

> Thiago -
> With regards to your original attempt for part b, I  believe using length()
> instead of sum() will get you what you want.  In this case, sum(x) and
> sum(x[x>1]) are the same assuming you have no values less than 1. Length can
> be used as a proxy for counting how many days rain>1.
>
>
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Count-values-greater-than-in-rasterbrick-layers-tp7586746p7586782.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

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