纬度和经度的栅格化

Question

给定raster对象r，如何创建一个新的raster，具有相同的范围和分辨率，并且单元格的值等于在r中相应单元格的纬度（或经度）？

例如，r可能如下所示：

r <- raster(matrix(runif(100), ncol=10))

- I Del Toro

---

你可以先创建/读取栅格，然后使用xFromCell()或yFromCell()提取每个单元格的纬度或经度，最后使用setValues()将提取的值分配为新的单元格值。 - shekeine

不要忘记标记语言 - 如果没有您的其他帖子，我们将不知道您正在尝试哪个环境。 （我已将其标记为 [tag:R]） - jbaums

我冒昧澄清了您的问题并添加了一个典型的栅格示例...希望您不介意。 - jbaums

Answer 1

您可以使用init函数。

library(raster)
r <- raster(nrow=10, ncol=10)
lon <- init(r, 'x')
lat <- init(r, 'y')

或者使用 terra

library(terra)
x <- rast(nrow=10, ncol=10)
lon <- init(x, 'x')
lat <- init(x, 'y')

- Robert Hijmans

---

整洁 - 我不知道这个。 - jbaums

Answer 2

一个简单的方法是(1) 复制栅格 r, (2) 用 coordinates提取其坐标，并将经度或纬度赋给新栅格对象的单元格。

例如，使用您的r:

library(raster)
r <- raster(matrix(runif(100), ncol=10))
lat <- lon <- r
xy <- coordinates(r)
lon[] <- xy[, 1]
lat[] <- xy[, 2]

以下是它们的外观:

plot(setNames(stack(r, lon, lat), c('r', 'lon', 'lat')))

- jbaums

---

太棒了！现在如果我想要整个地球的这张地图，我应该在哪里设置范围？ - I Del Toro

@IDelToro：尝试使用 r <- raster(extent(-180, 180, -90, 90), res=1)（在此情况下，将分辨率设置为以度为单位表示的首选单元格大小）。 - jbaums

Answer 3

如果您的问题是关于创建一个新的栅格对象，其范围和分辨率与另一个栅格对象相同，您可以使用命令template template是用于设置范围（以及在Raster*对象的情况下设置CRS）的Raster*或Extent对象。 如果不为NULL，则忽略参数xmn、xmx、ymn、ymx和crs（除非模板是Extent对象）。

r <- raster(matrix(runif(100), ncol=10))
r1 <- raster(x, template=r)

- Eko Susilo

Answer 4

如果有人只是想下载全球的纬度栅格图，并且不太关心它是如何制作的，我在这里发布了一个：

正负值（90S为-90）： https://drive.google.com/open?id=1fV1pnvlzi2PTJM7e6zx0S5IMP1mg09m-

仅正值（90S为90）： https://drive.google.com/open?id=1ibyNAp1c0E_DY9bFF-1gJw0vizKRmDUT

这些是以0.1度间隔发布的。

生成的Python代码：

import rasterio
import numpy as np

cellsize = .1
xedges = np.arange(-180,180,cellsize)
yedges = np.arange(90,-90,-cellsize)
XI, YI = np.meshgrid(xi,yi)

t = rasterio.transform.from_origin(xedges[0], yedges[0], cellsize, cellsize)

with rasterio.open('latitude.tif', 'w', driver='GTiff', 
                         height=YI.shape[0], width=YI.shape[1],
                         count=1, dtype=np.float32, transform=t) as src:
src.write(YI.astype(np.float32), 1)

with rasterio.open('latitude_abs.tif', 'w', driver='GTiff', 
                         height=YI.shape[0], width=YI.shape[1],
                         count=1, dtype=np.float32, transform=t) as src:
src.write(np.abs(YI).astype(np.float32), 1)

- Thomas Pingel