3

Sorry for cross-posting. I also posted this question at the R-SIG-GEO discussion list, but since I wanted to get as much feedback as possible I decided to post it here too.

I am trying to extract temperature values from a raster stack for about 400 municipalities in Brazil. My final goal is to create a data frame that is going to be used as a database for an interactive map server - probably using shiny and leaflet.

The final data frame would look like this:

> head(df)
Location        Var   Cut Year Month Freq
Campinas  temperature  10 2010  1    11
Campinas  temperature  10 2010  2    19
Campinas  temperature  10 2010  3    30
Campinas  temperature  10 2010  4    29
Campinas  temperature  10 2010  5    31
Campinas  temperature  10 2010  6    30

I have global raster stacks with daily data and I am counting, for each month in the raster, the number of days above certain temperature threshold. Please see below:

library(raster)
library(zoo)
library(maptools)

# Create a rasterStack similar to my data - same dimensions and layer names
r <- raster(ncol=360, nrow=180)
s <- stack(lapply(1:730, function(x) setValues(r, runif(ncell(r),min=0,max=30))))
idx <- seq(as.Date("2010/1/1"), by = "day", length.out = 730)
s <- setZ(s, idx)
s

# Define functions for 10, 15, 20 and 25 degrees
fun1 <- function(x, na.rm) {
sum(x > 10, na.rm)
}

fun2 <- function(x, na.rm) {
sum(x > 15, na.rm)
}

fun3 <- function(x, na.rm) {
sum(x > 20, na.rm)
}

fun4 <- function(x, na.rm) {
sum(x > 25, na.rm)
}

# Count number of days above the threshold temperature
days.above.10 <- zApply(s, by=as.yearmon, fun = fun1)
days.above.15 <- zApply(s, by=as.yearmon, fun = fun2)
days.above.20 <- zApply(s, by=as.yearmon, fun = fun3)
days.above.25 <- zApply(s, by=as.yearmon, fun = fun4)

Now, what I would like to do is to programmatically extract values for each location on my study area. The locations are defined as a shapefile with municipal contours of the Sao Paulo state in Brazil.

In this example, however, just for reproducibility's sake, I will be using a world polygon. But keep in mind that in my actual data the polygons will be much smaller.

# Import *sample* polygon data and subset only five "locations"
data(wrld_simpl)
locs <- subset(wrld_simpl, wrld_simpl@data$NAME %in% c("Argentina","Bolivia","Brazil","Paraguay","Uruguay"))

# Plot
plot(days.above.10,1)
plot(locs,add=T)

I feel like half of the work is done, but I am just grasping with the conversion to data frames.

Based on this self-contained example I provided, what would be the best strategy to come out with a data frame per location, like this?

> head(Argentina.df)
Location        Var Cut Year Month Freq
Argentina temperature  10 2010    1  11
Argentina temperature  10 2010    2  19
Argentina temperature  10 2010    3  30
Argentina temperature  10 2010    4  12
Argentina temperature  10 2010    5  17
Argentina temperature  10 2010    6  14

> head(Bolivia.df)
Location        Var  Cut Year Month Freq
Bolivia  temperature  10 2010    1  29
Bolivia  temperature  10 2010    2  31
Bolivia  temperature  10 2010    3  30
Bolivia  temperature  10 2010    4  17
Bolivia  temperature  10 2010    5  19
Bolivia  temperature  10 2010    6  12

and so on.

Note that "cut" refers to the temperature thresholds defined in the functions above. Each cut should come from the equivalent raster stack: days.above.10, days.above.15 and so on.

6

Since you are dealing with a rasterStack type of data you need to extract values using the extract function.The outcome therefore will be a vector.

As you see in the following line, the method for the extract function is:

extract(x,y,...)

in the position of x you will put the raster object and in the y the locations. Then use the as.data.frame function to turn the extracted object, into a dataframe.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.