在R中将矩阵转换为栅格图像

11

我有一个带有空间坐标和一个变量的矩阵数据。空间分辨率为1000米。

> 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 ...

我希望将其转换为光栅格式。

> 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

我使用以下代码来实现。但是我得到的分辨率与我实际数据的不同。有没有更好的方法来获得与我的真实数据相同的分辨率?

> 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)

如果您提供更小的示例代码,这些问题将更容易被人们调试。 - Spacedman
2个回答

21

试着阅读有关“光栅”的帮助。从矩阵创建光栅时,行和列的意义并不是您想象中的那样。您输入了一个大小为1241×710的矩阵,但从错误的向量中获取最大值和最小值。

请尝试以下操作:

> # 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)

好的,方形像素。现在创建你的光栅图:

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)

完全非正方形像素。如果你仔细看,矩阵从图像绘制中旋转了90度。或者说是转置之类的。

解决方法:只需从x、y、z列表中创建栅格:

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

像素是正方形的,与图像绘制方式一致,分辨率现在符合您的期望:

> 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)

10
既然我刚刚花了4个小时来写这个回答,那我再补充一下。当Spacedman说“行和列的概念不是你想象的那样”时,他的意思是,栅格图层是按行组织的,即从左到右,从上到下,而矩阵是按列组织的,即从上到下,从左到右。因此,在这种情况下,dat1[15]将给出与r[15]不同的结果。我个人认为栅格手册中并没有清晰地解释这一点。 - Dominik
1
你也可以转置矩阵(m): m <- t(m); raster(m[nrow(m):1,])。来源: https://gis.stackexchange.com/questions/60387/create-raster-file-in-r-from-matrix-and-latlons - Mikko

0

我认为你可能误解了这个上下文中“分辨率”的含义。它只是行和列中单元格数量的范围除以这个范围的大小。例如,在你的情况下,分辨率可以这样计算:

(1011340 - 302339.6) /1241
# 571.3138

(6671470  - 5431470) / 710
# 1746.479

所以,它的真正意思是每行有1241个正方形,由于每行从1011340302339.6单位(我认为在这种情况下是米),因此每个矩形的宽度为571。同样,每个矩形从上到下的距离为1746。

如果您希望每个光栅点都是边长为1000米的正方形,则横向范围不应该是(1011340-302339.6),而应该是1241*1000吗?


分辨率应该按照以下方式计算:'> (1011340-302339.6)/710 [1] 998.5921
(6671470 - 5431470)/1241 [1] 999.1942
- Jian Zhang

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接