全球地理距离栅格

Question

我在想是否有人已经建立了一张世界大陆的栅格图，其中每个单元格等于该单元格到最近海岸的距离。这张地图将突出显示最为孤立的内陆地区。

我想这只需要对全球边界的shapefile进行光栅化，然后计算出距离即可。

- I Del Toro

Answer 1

您可以使用raster::distance函数来实现此操作，该函数计算每个NA单元格到最近的非NA单元格的距离。您只需要创建一个栅格图像，其中陆地像素为NA，而非陆地像素则为其他值。

具体方法如下：

library(raster)
library(maptools)
data(wrld_simpl)

# Create a raster template for rasterizing the polys. 
# (set the desired grid resolution with res)
r <- raster(xmn=-180, xmx=180, ymn=-90, ymx=90, res=1)

# Rasterize and set land pixels to NA
r2 <- rasterize(wrld_simpl, r, 1)
r3 <- mask(is.na(r2), r2, maskvalue=1, updatevalue=NA)

# Calculate distance to nearest non-NA pixel
d <- distance(r3)

# Optionally set non-land pixels to NA (otherwise values are "distance to non-land")
d <- d*r2

为了创建上面的图形（我喜欢使用rasterVis进行绘图，但您也可以使用plot(r)）：

library(rasterVis)
levelplot(d/1000, margin=FALSE, at=seq(0, maxValue(d)/1000, length=100),
          colorkey=list(height=0.6), main='Distance to coast')

- jbaums

---

很棒的回答！两个后续问题：1）像素的分辨率可以进一步降低吗？2）地图的值被解释为公里？ - I Del Toro

@IDelToro（1）您可以设置所需的分辨率，其中我们定义“r”（此处为度数中的“res”参数）; （2）查看?raster-如果数据是“未投影的”（即地理空间），则distance的输出以米为单位。否则，它以投影单位为单位。请注意，在绘制时，我通过除以1000来得到公里。 - jbaums

谢谢！您会用相同的方式创建一个纬度（或经度）值的栅格吗？ - I Del Toro

@IDelToro：你可以从其他栅格数据中派生这些，例如 lon <- lat <- d; lon[] <- coordinates(d)[, 1]; lat[] <- coordinates(d)[, 2]。如果要从头开始创建，请使用矩阵 m，例如经度，然后使用 raster(m)。 - jbaums

感谢您的回答。一般来说，我发现将分辨率增加到0.5会显著增加处理时间（这很好），但进一步增加到0.1度时永远无法完成。还有其他人遇到过这种情况吗？没有错误提示，它只是占用了一个核心的100％，并没有提示它是否正在运行或已经挂起。 - dez93_2000