3

I have a bunch of raster files in different formats: - NETCDF 4 - HDF 5 - GeoTIFF

Using rgal package I'm able to open any of them and use them as objects in R! Creating stacks of them allows me to calculate correlation coefficient.

But what if I would like to mąkę some calculations on one pixel only (for example to create scatterplot)?

Am I able to extract value for particular coordinates? Obviously I can use gdal2xyz tool, but I would like to do that entirely in R!

4
  • If you do not mind trying another package, {raster} package provides getValues() function to extract cell values of a specified row. You can then narrow down column by [ ]. For example getValues(your_raster_object, 10)[20] will give you the value of cell at (row= 10, column= 20). – Kazuhito Nov 13 '16 at 11:51
  • 1
    Use raster::brick to open, then extract() to get time series at point locations (or cell indexes) – mdsumner Nov 13 '16 at 12:34
  • @Kazuhito I think that you're reffering to single raster, not raster stack object. Anyway, could you post your comments as responses, I think that they're useful. – matandked Nov 19 '16 at 14:22
  • @matandked Thanks pointing out that, sorry. Will post it as-is; please add your comment on it. Cheers, – Kazuhito Nov 19 '16 at 22:49
2

You can use as.data.frame(your_raster_data, xy = T) to convert you data into a data.frame and plot it. If layers represents different dates, convert your data to zoo object:

library(raster)
library(zoo)
r <- raster()
r <- setValues(r, values = 1:length(r))
r2 <- r*2
r3 <- r*4
r4 <- r*3

r.stack <- stack(r,r2,r3,r4)

names(r.stack) <- c("D_2016_11_01","D_2016_11_06","D_2016_11_07","D_2016_11_10")

df <- as.data.frame(r.stack, xy = T)

date <- as.Date(substr(names(r.stack),3,12), format = "%Y_%m_%d")

to.plot <- zoo(unlist(df[df$x == 45.5 & df$y ==85.5,c(3:dim(df)[2])]), date)

plot(to.plot)

enter image description here

2
  • @matandked To filter by value in a data.frame, you need to use data.frame[row condition,column condition]. Remember this. – aldo_tapia Nov 20 '16 at 10:46
  • Yes, sorry for this very stupid mistake - I've tried different approaches, even without comma which is obviously incorrect I finally manage how to accomplish what I want - create scatterplot for selected pixel: gist.github.com/kedziorm/a43d413d7cedd004815171ed32db9db8 – matandked Nov 26 '16 at 22:47
1

If you do not mind trying another package, {raster} package provides getValues() function to extract cell values of a specified row. You can then narrow down column by [ ]. For example getValues(your_raster_object, 10)[20] will give you the value of cell at (row= 10, column= 20).

[EDIT]

The question was about accessing RasterBrick cells. Please try following example and find how [[ ]] and [ ] would work with getValues().

filename <- system.file("external/rlogo.grd", package="raster")
(b <- brick(filename))
b[[1]]    # NB. you will obtain RasterLayer
b[[2:3]]  # NB. you will obtain RasterStack

getValues(b[[1]], 77)       # Layer:1, Row 77
getValues(b[[1]], 77)[101]  # Layer:1, Row 77, Column 101
5
  • I have following RasterBrick: class : RasterBrick dimensions : 19, 33, 627, 55 (nrow, ncol, ncell, nlayers) resolution : 0.259366, 0.3138765 (x, y) extent : 16.85879, 25.41787, 48.69525, 54.65891 (xmin, xmax, ymin, ymax) When I type 'getValues(myRasterBrick,18)' I get numeruous values for all rasters, but when I type 'getValues(SMOS,18)[1]' I get only one value. Moreover, how I will now which row/column corresponds to which coordinate value? – matandked Nov 20 '16 at 7:27
  • @matandked Hi, I see multiple questions which would fit for other posts... If I can comment just about accessing individual layer, I think you can use double square brackets [[ to get subset of bands. (e.g. [[2]] or [[1:3]]). – Kazuhito Nov 20 '16 at 8:22
  • @matandked My apologies... you were clarifying the point that you were making in your original question for me... I will revise my answer to reflect it. As to your last question, it will make a good new post. – Kazuhito Nov 20 '16 at 9:26
  • Thank you for your responses, however, I still have an issue with retrieving value. I've tried to follow your example, but when I type: nlayers = length(myRasterBrick[1]) getValues(myRasterBrick[[1:nlayers]],1)[1] I receive only one value - I would like to select one pixel and retrieve values for this particular pixel from all raster layers. In other words, I would like to retrieve the vector of values. Length of the vector should be equal to the number of layers. – matandked Nov 26 '16 at 21:15
  • @matandked Understood you are interested in obtaining cell values across the layers of the brick. Then getValues() is not the right one. Using the same example, try getValuesBlock(b, row=35, nrows=3, col=50, ncols=3, lyrs=1:3). If you need not to control which layers to choose, it is even easier. Then try b[35, 50][1, ]. Cheers, – Kazuhito Nov 27 '16 at 1:16

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.