How to read and write raster matrix data with R (raster package and image processing)?

Question

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.

Answer 1

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)

Answer 2

Not exactly what you asked for, but this whole raster-to-matrix and matrix-to-raster conversion can also be performed making use of {terra}, slowly replacing {raster}:

library(terra)
#> terra 1.6.17

r1 <- system.file("ex/elev.tif", package="terra") |> rast()
r1
#> class       : SpatRaster 
#> dimensions  : 90, 95, 1  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : elev.tif 
#> name        : elevation 
#> min value   :       141 
#> max value   :       547

m <- as.matrix(r1, wide = TRUE)
dim(m)
#> [1] 90 95

# --- some black box operations on your matrix object ---

r2 <- rast(m, crs = "epsg:4326", ext = ext(r1))
r2
#> class       : SpatRaster 
#> dimensions  : 90, 95, 1  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : memory 
#> name        : lyr.1 
#> min value   :   141 
#> max value   :   547

^{Created on 2022-10-14 with reprex v2.0.2}

As you can see, extent, resolution and coordinate reference system are preserved, so that no unintentional shifting is to be expected here.

In order to write the SpatRast object to disk, you can make use of terra::writeRaster().