Read and write raster matrix data with R (raster package and image processing)

digital image processingmatrixrraster

I'm using R for morphological operations on geo-rasters.
The packages mmnand and EBImage offer morphological operations. However they can't work with georeferenced rasters, but operate on a matrix (or their own classes).

So in order to use them, I need to (1) "extract" the raster data as matrix, (2) manipulate it with EBImage (being slightly faster than mmnand) and then (3) store it back to the raster, to reestablish the geospatial information (CRS, coordinates, …) and writeRaster as GeoTiff.

This works, but somehow messes with the coordinates/positioning. The whole raster is shifted see here for examople image of the shift (Imgur image upload currently broken).

How can I correctly read raster data as matrix and reassign it to a raster?

Here is code for what I'm doing currently:

# read file
mask_raster = raster("my_file.tif")

# (1) "extract" matrix
mask_matrix = as.matrix(mask_raster) 

# (2) make kernel and perform morphological operations
kernel = makeBrush(3, shape='box')
mask_matrix = opening(mask_matrix, kernel)
mask_matrix = closing(mask_matrix, kernel)

# (3) store matrix data back to raster to "restore" geospatial context 
mask_raster[] = mask_matrix

# write geotiff
writeRaster(mask_raster, "my_file_morphed.tif", datatype ="INT1U", options="COMPRESS=LZW", overwrite = T)

I presume the key/problem is, that as.matrix (an array?) gives something different then my_raster[] (a vector?), but I'm not familiar enough with the raw data structures to determine how to solve this "incompatibility".

  mask_matrix = as.matrix(mask_raster) 
  # creates a "Large matrix (2903145 elements, 11.1 Mb): int [1:1505, 1:1929]"

  mask_vector = mask_raster[] 
  # creates a "Large integer (2903145 elements, 11.1 Mb): int [1:2903145]"

To perform the morphological operations, I need the matrix. When assigning it via mask_raster[] = mask_matrix I tried to convert it to vector mask_raster[] = c(t(mask_matrix)) (needs transposing with t(), otherwise the row/col order is swapped), but this doesn't solve the shift. Comparing the two raster objects (before and after assignment of the "new" vector) doesn't show any differences… I'm really puzzled here.

Best Answer

I think what you want to do is use m = as.matrix(r) to get a matrix from a raster, and r2 = raster(m);extent(r2) = extent(r) to get a raster from a matrix, which you then have to set the extent on to get a compatible raster.

Here's a worked example that shows a round trip from matrix to raster to matrix to raster:

First let's make a raster from this matrix:

> m = matrix(1:12, 3, 4)
> r = raster(m)

and we'll give it an extent:

> extent(r)=c(1,2,3,4)
> r
class      : RasterLayer 
dimensions : 3, 4, 12  (nrow, ncol, ncell)
resolution : 0.25, 0.3333333  (x, y)
extent     : 1, 2, 3, 4  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : 1, 12  (min, max)

Now we get the matrix of that, do something with it (that importantly preserves the number of rows and columns) that is like whatever image processing you are doing in the matrix world:

> mm = as.matrix(r)
> mm
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
> mm2 = mm*2

Now to get mm2 back in the raster world, convert from matrix:

> rmm2 = raster(mm2)
> rmm2
class      : RasterLayer 
dimensions : 3, 4, 12  (nrow, ncol, ncell)
resolution : 0.25, 0.3333333  (x, y)
extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : 2, 24  (min, max)

and then set the extent:

> extent(rmm2) = extent(r)
> rmm2
class      : RasterLayer 
dimensions : 3, 4, 12  (nrow, ncol, ncell)
resolution : 0.25, 0.3333333  (x, y)
extent     : 1, 2, 3, 4  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : 2, 24  (min, max)

You can do this conversion in one step but its a lot of typing so put it in a function so you only have to do it once:

> rmm2 = raster(mm2, xmn=xmin(r),xmx=xmax(r),ymn=ymin(r), ymx=ymax(r))
> rmm2
class      : RasterLayer 
dimensions : 3, 4, 12  (nrow, ncol, ncell)
resolution : 0.25, 0.3333333  (x, y)
extent     : 1, 2, 3, 4  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : 2, 24  (min, max)
Related Question