convert matrix to raster in R
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)
Jian Zhang
Updated on July 12, 2022Comments
-
Jian Zhang almost 2 years
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)
-
Jian Zhang over 11 yearsThe resolution should be calculated like this: '> (1011340-302339.6)/710 [1] 998.5921 > (6671470 - 5431470)/1241 [1] 999.1942
-
Dominik over 9 yearsSince 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.
-
Mikko over 4 yearsYou can also transpose the matrix (m):
m <- t(m); raster(m[nrow(m):1,])
. Source: gis.stackexchange.com/questions/60387/…