extract raster pixel value as well as pixel xy coordinates using a SpatialLine object in R

Question

I have a Spatialline that I converted from a polygon shapefile (digitized manually based on features in "imagebrick" - this means that spatially the "polyline" and "imagebrick" are overlapped as I wanted)

polyline <- as(shapefiles_data[1,],"SpatialLines")

> polyline
class       : SpatialLines 
features    : 1 
extent      : 357714.3, 357719, 4076030, 4076035  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=11 +datum=NAD27 +units=m +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat

And a Rasterbrick

> imagebrick
class       : RasterBrick 
dimensions  : 29180, 14940, 435949200, 4  (nrow, ncol, ncell, nlayers)
resolution  : 0.6, 0.6  (x, y)
extent      : 354038.4, 363002.4, 4066992, 4084500  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=11 +datum=NAD27 +units=m +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat 
names       : band1, band2, band3, band4 
min values  :     0,     0,     0,     0 
max values  : 65535, 65535, 65535, 65535

I was trying to extract pixels in the Rasterbrick that touched by the Line(SpatialLine). But I also wanted to pull out xy coordinates associated with those extracted pixels and store them together in a SpatialPointsDataFrame. I've tried:

extract_polyline_test <- extract(imagebrick, polyline)
extract_polyline_test <- as.data.frame(extract_polyline_test)
> extract_polyline_test
band1 band2 band3 band4
1    125   156    73   392
2    129   161    80   366
3    156   208   122   374
4    142   175    97   330
5    102   117    31   389
6    117   143    57   381
7    162   221   127   413
8    213   302   204   405
..   ...   ...   ...   ...

Then, I was trying to use RasterToPoints function to pull out the xy coordinates. This requires a raster layer that contains all the pixels that I want xy from. I've tried to use "crop" to get that raster "line" but it does not work.

rastercrop_polyline <-crop(imagebrick,polyline)
rastercrop_polyline

> rastercrop_polyline
class       : RasterBrick 
dimensions  : 8, 7, 56, 4  (nrow, ncol, ncell, nlayers)
resolution  : 0.6, 0.6  (x, y)
extent      : 357714.6, 357718.8, 4076030, 4076035  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=11 +datum=NAD27 +units=m +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat 
data source : in memory
names       : band1, band2, band3, band4 
min values  :    72,    63,     0,   156 
max values  :   274,   415,   299,   497 

s.polyline <- rasterToPoints(rastercrop_polyline)
> s.polyline
         x       y band1 band2 band3 band4
 [1,] 357714.9 4076035    93    96    17   342
 [2,] 357715.5 4076035    72    63     0   377
 [3,] 357716.1 4076035    90    93    17   371
 [4,] 357716.7 4076035   125   156    73   392
 [5,] 357717.3 4076035   129   161    80   366
 [6,] 357717.9 4076035   156   208   122   374
 [7,] 357718.5 4076035   142   175    97   330
 [8,] 357714.9 4076034   100   108    25   326
 ..  ........  .......   ...   ..     ..    ...

It seems that the "crop" function does not crop a "line" of raster but a whole "polygon". Is there any other way to get the xy coordiantes of the pixels that I extracted?

p.s. I was also trying to extract surrounding pixels of each extracted pixels. How should I do this? I notice that there is a function parameter in extract function. Should I work on writing a small function within the extraction function or there are any other ways to do this?

p.p.s I was also thinking to use Rasterize function to rasterize the "SpatialLine" and then use the RasterToPoints to pull out the xy coordinates. But I can not control each pixel of the "rasterized SpatialLine" equals to the pixels that touched by the "line" in the original imagebrick.

Thanks everyone for your time viewing my question and thanks in advance for your help.

Robert Hijmans Robert Hijmans · Accepted Answer · 2014-03-18T04:48:03

You say that you were "thinking to use Rasterize function to rasterize the "SpatialLine" and then use the RasterToPoints to pull out the xy coordinates."

That seems correct for what you want

But then you say that you "can not control each pixel of the "rasterized SpatialLine" equals to the pixels that touched by the "line" in the original imagebrick."

But they are, except that cells with multiple lines are only represented once.

Another approach:

# set up some example data
r <- raster()
values(r) <- runif(ncell(r))
cds1 <- rbind(c(-50,0), c(0,60), c(40,5), c(15,-45), c(-10,-25))
cds2 <- rbind(c(80,20), c(140,60), c(160,0), c(140,-55))
lines <- SpatialLines(list(Lines(list(Line(cds1)), "1"), Lines(list(Line(cds2)), "2") ))

# values for lines
v <- extract(r, lines)
# create raster with cell numbers
ids <- init(r, v='cell')
# extract these
cells <- extract(ids, lines)
# compute xy
xy = lapply(cells, function(x) xyFromCell(r, x))

extract raster pixel value as well as pixel xy coordinates using a SpatialLine object in R

1 Answers