Extract values along customized circular strips (constant radial distance) from raster

Question

I have a large raster object which was imported from a .tif image file. I want to analyze pixels of this raster, at a given radial distance from the center, in order to identify certain axisymmetric phenomenon that I can notice in the image (below). To do this I want to extract the values of pixels that intersect with circular strips of a given radii (and width) from the center of the image (indicated in the figure below).

I have explored several options to do this including the extract function and the imager package. In the imager package, you can easily extract values along rows or columns but I couldn't find functions to extract values with customized shapes like I require.

I might use the extract function with SpatialPolygons argument, however for this I will have to input the locations of all the points in a polygon object, and this would require very high density location of points (since the points, I think, are joined by line elements). Moreover, I want to vary the number density of the strips over which I extract pixel values (and later average) therefore this method would be both tedious and inflexible. I was thus wondering if any of you have any suggestions to tackle this issue.

> str(imRaster)
Formal class 'RasterLayer' [package "raster"] with 12 slots
  ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
  .. .. ..@ name        : chr "C:\\Users\\Nandu\\input_images\\All\\101_1A_1000ms.tif"
  .. .. ..@ datanotation: chr "INT2U"
  .. .. ..@ byteorder   : chr "little"
  .. .. ..@ nodatavalue : num -Inf
  .. .. ..@ NAchanged   : logi FALSE
  .. .. ..@ nbands      : int 1
  .. .. ..@ blockrows   : int 1
  .. .. ..@ blockcols   : int 1392
  .. .. ..@ driver      : chr "gdal"
  .. .. ..@ open        : logi FALSE
  ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
       .
       .
       .

Alternately, if there are other methods (other than by importing a .tif image to raster) that would allow such an operation, that might also be useful. Please let me know.

Thanks in advance!

EDIT: In order to make the problem reproducible, I have further added here the link to the image that I want to use. It is a tiff image, and can be downloaded from here

Link to the tiff image that I am trying to process in R

library(tiff)
library(raster)
imRaster = raster(file_path to the image)
plot(imRaster)
xy = cbind(684.4228, 599.0458) # Gives a rough location of the center in the image coordinates

This code can be run in order to obtain the raster described above and visualize it. Hope it helps!

This is an interesting question. Please provide a [minimal, reproducible example](http://stackoverflow.com/help/mcve). Therefore, you might have a look at [how to write a good reproducible R example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) — loki, May 04 '17 at 13:00
Thanks @loki. I have added a link to the tiff image and a sample code to get started immediately. Hope it helps. Let me know if anything else is needed. — Pradical2190, May 05 '17 at 10:26

loki · Accepted Answer · 2017-05-07T17:51:17.910

Let's start by loading the raster and assigning it a fake coordinate reference system (CRS). We use a metrical one so that 1 pixel in your image refers to one "meter" in the pseudo-geographic raster.

library(raster)

imRaster <- raster("~/Downloads/46-ECN300448-216.tif")
imRaster@crs <-  CRS("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs")

plot(imRaster, col = grey.colors(255))

Afterwards, as you mentioned we use your center point xy and convert it to a SpatialPoint

xy <- cbind(684.4228, 599.0458)
xySp <- SpatialPoints(coords = xy)
plot(xySp, add = T)

Then we use the rgeos package to create two radial buffers around the point. In this example we use 200 and 300 px/m.

library(rgeos)

outer <- 300
inner <- 200

bufOut <- gBuffer(xySp, width = outer)
bufIn <- gBuffer(xySp, width = inner)

strip <- bufOut - bufIn
plot(strip, add = T, col = "#FF000050")

Finally, we can use the strip to mask the raster image and calculate statistics (or even use the original values) with getValues().

m <- mask(imRaster, mask = strip)

# plot(m) # plot the mask if you want to see what it looks like

mean(getValues(m), na.rm = T)
# [1] 17004.24

Important Note: Only the cells having their center within the polygon are remain after masking. You might tackle this issue by playing around with the with of the buffers.

Nice approach. However I can't generate the strip simply by subtracting the two geometries. I would have used the function `gDifference` to obtain the strip object... Is it just me? — G. Cocca, May 06 '17 at 14:25
@G.Cocca, I checked this approach both on a Windows and an Ubuntu (16.04 LTS) System. Newest R and package versions. Have you tried updating your packages? However, `gDifference` also does the trick. @Pradical2190, does this solution satisfy you? — loki, May 07 '17 at 17:54
@loki, thanks for the solution. It works well! I also wanted to note, before this answer, I was hunting for ways to do this and realized that this can be achieved by using the extract function from the raster package itself, and using the buffer argument. However, your method has the advantage that you get the selected regions as a SpatialPolygons class, and that is more convenient to manipulate later. Hence I have also accepted your answer. Thanks! — Pradical2190, May 09 '17 at 09:25

Extract values along customized circular strips (constant radial distance) from raster

1 Answers1