How to Extract Data from a Rasterbrick

How to extract data from a RasterBrick?

lets take a grid of 3x4 rasters over three years in a silly calendar that only has seven months in it:

d = array(1:(3*4*7*3),c(3,4,7*3))
b = brick(d)

Now lets give the brick layers names by year and month:

names(b) = paste("rain",outer(1:7,2001:2003,paste,sep="-"),sep="-")
> names(b)
 [1] "rain.1.2001" "rain.2.2001" "rain.3.2001" "rain.4.2001" "rain.5.2001"
 [6] "rain.6.2001" "rain.7.2001" "rain.1.2002" "rain.2.2002" "rain.3.2002"
[11] "rain.4.2002" "rain.5.2002" "rain.6.2002" "rain.7.2002" "rain.1.2003"
[16] "rain.2.2003" "rain.3.2003" "rain.4.2003" "rain.5.2003" "rain.6.2003"
[21] "rain.7.2003"

and make some test points:

> pts = data.frame(x=runif(3),y=runif(3), month=c(5,1,3),year = c(2001,2001,2003))
> pts
          x         y month year
1 0.2513102 0.8552493     5 2001
2 0.4268405 0.3261680     1 2001
3 0.7228359 0.7607707     3 2003

Now construct the layer name for the points, and match to the names:

pts$layername = paste("rain",pts$month,pts$year,sep=".")
pts$layerindex = match(pts$layername, names(b))

Now I don't think the layer index in extract is vectorised, so you have to do it in a loop...

> lapply(1:nrow(pts), function(i){extract(b, cbind(pts$x[i],pts$y[i]), layer=pts$layerindex[i], nl=1)})
[[1]]
     rain.5.2001
[1,]          57

[[2]]
     rain.1.2001
[1,]           5

[[3]]
     rain.3.2003
[1,]         201

Or in a simple vector:

> sapply(1:nrow(pts), function(i){extract(b, cbind(pts$x[i],pts$y[i]), layer=pts$layerindex[i], nl=1)})
[1]  57   5 201

I'd do some checks to make sure those values are what you expect from those inputs before doing it on anything major though. Its easy to get indexes the wrong way round....

Another way to do it with a single extract call is to compute the values for all layers and then extract with a 2-column matrix subset:

> extract(b, cbind(pts$x, pts$y))[
      cbind(1:nrow(pts),match(pts$layername, names(b)))
     ]
[1]  57   5 201

Same numbers, comfortingly.

Extracting data from lower layers in a Rasterbrick

Your code is taking the value from the layer of the previous point, not the previous layer.

To see that imagine we are looking at the point in row 2 (i=2). your code that indicates the layer is pts$layerindex[i-1], which is pts$layerindex[1]. In other words, the layer of the point in row 1.

The fix is easy enough. For clarity I will write the function separetely:

foo = function(i) extract(b, cbind(pts$x[i],pts$y[i]), layer=pts$layerindex[i]-1, nl=1)
sapply(1:nrow(pts), foo)

I have not tested it, but this should be all.

Extracting data for certain coordinates from RasterBrick

How about using tidyr::gather?

xy <- raster::extract(rr, cbind(9.258157, 47.70842), df = TRUE)[, 1:10]

library(tidyr)
xy
  ID X1950.01.01 X1950.01.02 X1950.01.03 X1950.01.04 X1950.01.05 X1950.01.06 X1950.01.07 X1950.01.08 X1950.01.09
1  1           0         9.3         4.2        11.2           0           0           0           0           0

gather(xy[, -1], key = datum, value = Prcp)
        datum Prcp
1 X1950.01.01  0.0
2 X1950.01.02  9.3
3 X1950.01.03  4.2
4 X1950.01.04 11.2
5 X1950.01.05  0.0
6 X1950.01.06  0.0
7 X1950.01.07  0.0
8 X1950.01.08  0.0
9 X1950.01.09  0.0

Extracting all individual layers from a Raster Brick File

writeRaster seems to strip off the dot-number before creating a raster file. Hence it tries to write your layers all to Data.Fiel//tNO2Trop.tif.

> writeRaster(r, "./test.2", format="GTiff")
> dir(".")
[1] "test.tif"

(Note for some reason your code has format("GTiff") for format="GTiff". This works by the fluke that format is a function and returns the string "GTiff" and writeRaster is expecting the format string here)

I don't know why and I don't know if this is documented or a bug. You can work round by using dashes instead of dots:

> writeRaster(r, "./test-2", format="GTiff")
> dir(".")
[1] "test-2.tif" "test.tif" 

and if dots are important to you then do a file.rename afterwards.

Edit: If you add the .tif to the file names then all is well:

> writeRaster(s, names(s), bylayer=TRUE, format="GTiff")
Error in .getGDALtransient(x, filename = filename, options = options,  : 
  filename exists; use overwrite=TRUE

fails on the second layer because dot-number is stripped:

> dir()
[1] "layer.tif"

add .tif to the names:

> writeRaster(s, paste0(names(s),".tif"), bylayer=TRUE, format="GTiff")

shazam:

> dir()
[1] "layer.1.tif" "layer.2.tif" "layer.3.tif" "layer.tif"  

Conditionally extract data from a raster stack data based on values in a spatial points dataframe

Here's one option:

library(raster)

# Vector of dates
dates <- format(seq(as.Date('1999/1/11'), as.Date('2006/1/10'), by='month'), '%Y%m')

# RasterStack with random data
s <- setNames(stack(replicate(length(dates), raster(matrix(runif(100), 10)))), 
              paste0('rain', dates))

# Create a SpatialPointsDataFrame with some random dates and coords
d <- data.frame(x=runif(10), y=runif(10), date=sample(dates, 10))
coordinates(d) <- ~x+y

# Split the spdf by date
d_by_date <- split(d, d$date)

# Extract values
rain_sum <- unsplit(lapply(d_by_date, function(x) {
  # current year
  y <- as.numeric(substr(x$date, 1, 4))
  # current month
  m <- as.numeric(substr(x$date, 5, 6))
  # if month is after Oct, start from that year's Oct
  # if month is before Oct, start from previous year's Oct
  if(m < 11) y <- y-1
  start_date <- as.Date(sprintf('%s/10/01', y))
  # if start_date is earlier than first time slice, reset to first time slice
  start_date <- max(min(as.Date(sub('rain', '01', names(s)), '%d%Y%m')), start_date)
  end_date <- as.Date(paste0(x$date, '01'), '%Y%m%d')
  # Sequence of dates to sum over
  i <- format(seq(start_date, end_date, by='month'), 'rain%Y%m')
  # Extract values and sum
  sum(extract(s[[i]], x))
}), d$date)

---

Related Topics