convert matrix to raster in R

Question

I have a matrix data with spatial coordinates and one variable. The spatial resolution is 1000 meters.

> str(dat1)
> List of 3
> $ x: num [1:710] 302340 303340 304340 305340 306340 ...
> $ y: num [1:1241] 5431470 5432470 5433470 5434470 5435470 ...
> $ z: num [1:710, 1:1241] 225 225 225 225 225 ...

I want to convert it into raster format.

> dat1$x[1:10]
> [1] 302339.6 303339.6 304339.6 305339.6 306339.6 307339.6 308339.6 309339.6 310339.6 311339.6
> dat1$y[1:10]
>  [1] 5431470 5432470 5433470 5434470 5435470 5436470 5437470 5438470 5439470 5440470

I used the following code to do it. But the resolution I get is not the same with that I have. Any better way to get the same resolution with my real data?

> r <-raster(
             dat1$z,
             xmn=range(dat1$x)[1], xmx=range(dat1$x)[2],
             ymn=range(dat1$y)[1], ymx=range(dat1$y)[2], 
             crs=CRS("+proj=utm +zone=11 +datum=NAD83")
            )
> r

class       : RasterLayer 
dimensions  : 710, 1241, 881110  (nrow, ncol, ncell)
resolution  : 571.3135, 1746.479  (x, y)
extent      : 302339.6, 1011340, 5431470, 6671470  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=11 +datum=NAD83 
data source : in memory
names       : layer 
values      : 13.65059, 248.6229  (min, max)

Aldo

These things are a lot easier for people to debug if you make smaller examples and give us the code to make them. – Nancee

Answer 1

Try reading the help for raster. When creating a raster from a matrix, the sense of rows and columns isn't what you think it is. You were feeding it a 1241x710 matrix but taking the max and min from the wrong vectors.

Try the following:

> # small version of your test set
> dat1=list()
> dat1$x=seq(302339.6,by=1000,len=71)
> dat1$y=seq(5431470,by=1000,len=124)
> dat1$z=matrix(runif(71*124),71,124)
> str(dat1)
List of 3
 $ x: num [1:71] 302340 303340 304340 305340 306340 ...
 $ y: num [1:124] 5431470 5432470 5433470 5434470 5435470 ...
 $ z: num [1:71, 1:124] 0.765 0.79 0.185 0.461 0.421 ...
> image(dat1,asp=1)

Nice square pixels. Now create your raster:

r <-raster(
             dat1$z,
             xmn=range(dat1$x)[1], xmx=range(dat1$x)[2],
             ymn=range(dat1$y)[1], ymx=range(dat1$y)[2], 
             crs=CRS("+proj=utm +zone=11 +datum=NAD83")
            )
plot(r)

Totally NON-square pixels. And if you look carefully, the matrix is rotated 90 degrees from the image plot. Or transposed or something.

Solution: just create the raster from the x,y,z list:

 > r=raster(dat1);plot(r)

Square pixels, same way round as image plot, and resolution is now what you expect:

> r
class       : RasterLayer 
dimensions  : 124, 71, 8804  (nrow, ncol, ncell)
resolution  : 1000, 1000  (x, y)
extent      : 301839.6, 372839.6, 5430970, 5554970  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : in memory
names       : layer 
values      : 7.738103e-05, 0.9995497  (min, max)

Nancee

Since I just spent 4 hours on this, I'll add to this answer. What Spacedman means when he says "the sense of rows and columns isn't what you think it is" he means that a rasterlayer is organized by row, i.e. left-right, top-bottom, while a matrix is organized by column, i.e. top-bottom, left-right. so in this case, dat1[15] will give you a different result than r[15]. I personally don't think this is spelled out clearly enough in the raster manual. – Perrin

You can also transpose the matrix (m): m <- t(m); raster(m[nrow(m):1,]). Source: gis.stackexchange.com/questions/60387/… – Dalhousie

Answer 2

I think you might misunderstand what resolution means in this context. It is simply the extent divided by the number of cells within the rows and columns. For example, in your case, the resolution was calculated like this:

(1011340 - 302339.6) /1241
# 571.3138

(6671470  - 5431470) / 710
# 1746.479

So, what it really means is that you have 1241 squares in each row, and since each row goes from 1011340 to 302339.6 units (I think metres in this case), then each rectangle is 571 across. Similarly, each rectangle is 1746 from top to bottom.

If you expect that each raster point to be a square with sides equal to 1000 metres, then shouldn't the range across be 1241*1000 (instead of (1011340 - 302339.6))?

Hoeg

The resolution should be calculated like this: '> (1011340-302339.6)/710 [1] 998.5921 > (6671470 - 5431470)/1241 [1] 999.1942 – Aldo