在不规则网格上绘制气候数据的正确方法

Posted 2023-02-16

【中文标题】在不规则网格上绘制气候数据的正确方法

【英文标题】：A proper way to plot climate data on an irregular grid

【发布时间】：2018-08-09 09:09:36

【问题描述】：

我已将这个问题作为Efficient way to plot data on an irregular grid 问题的一部分提出，但一般反馈是将原始问题拆分为更易于管理的部分。因此，这个新问题。

我使用在不规则二维网格上组织的卫星数据，其维度是扫描线（沿轨道维度，即 Y 轴）和地面像素（跨轨道维度，即 X 轴）。每个中心像素的经纬度信息存储在辅助坐标变量中，以及四个角坐标对（经纬度坐标在 WGS84 参考椭球上给出）。

让我们构建一个玩具数据集，其中包含一个 12x10 的潜在不规则网格和相关的表面温度测量值。

library(pracma) # for the meshgrid function
library(ggplot2)

num_sl <- 12 # number of scanlines
num_gp <- 10 # number of ground pixels
l <- meshgrid(seq(from=-20, to=20, length.out = num_gp), 
              seq(from=30, to=60, length.out = num_sl))

lon <- l[[1]] + l[[2]]/10
lat <- l[[2]] + l[[1]]/10

data <- matrix(seq(from = 30, to = 0, length.out = num_sl*num_gp), 
               byrow = TRUE, nrow = num_sl, ncol = num_gp) +
  matrix(runif(num_gp*num_sl)*6, nrow = num_sl, ncol = num_gp)

df <- data.frame(lat=as.vector(lat), lon=as.vector(lon), temp=as.vector(data))

lon 和 lat 数据包含我正在使用的原始产品中提供的中心像素坐标，存储为二维矩阵，其轴为 ground_pixel（X 轴）和扫描线（Y 轴） . data 矩阵（相同的维度）包含我的测量值。然后我展平这三个矩阵并将它们存储在一个数据框中。

我想在地图上绘制地面像素（作为四边形），并相应地填充温度测量值。

使用我得到的图块：

ggplot(df, aes(y=lat, x=lon, fill=temp)) + 
  geom_tile(width=2, height=2) +
  geom_point(size=.1) +
  borders('world', colour='gray50', size=.2) + 
  coord_quickmap(xlim=range(lon), ylim=range(lat)) +
  scale_fill_distiller(palette='Spectral') +
  theme_minimal()

但这不是我所追求的。我可以使用width 和height 来使图块相互“接触”，但这当然不会接近我想要的目标，即绘制实际的投影地图上的地面像素。 例如，Python 的 xarray 可以根据像素中心坐标自动推断像素边界：

问题

有没有办法在 R 中实现相同的结果，即：从像素中心自动推断像素边界，并将像素绘制为地图上的填充多边形？也许使用sf 包？

我可以在这个question 的答案中看到它已经完成，但是我认为使用sf 的答案有点不清楚，因为它处理不同的投影和潜在的规则网格，而在我的情况下，我想我不必重新投影我的数据，而且，我的数据不在常规网格上。

如果这是不可能的，我想我可以在我的产品中使用像素边界信息，但如果这个问题被证明不容易解决，那么这可能是另一个问题的主题。

【问题讨论】：

你说你有每个瓦片角的坐标？我建议使用sf 创建平铺网格并使用ggplot 的开发版本使用geom_sf 进行绘图。如果在制作这些多边形时正确设置了 CRS，则应该可以获得所需的 python 图。具体如何做到这一点取决于坐标和温度测量值的存储方式——当前示例数据只有中心像素，对吧？

是的，没错。我希望有一种简单的方法可以从像素中心推断像素边界。我已经看到了这个：polys = as(SpatialPixelsDataFrame(orig_grid, orig_grid@data, tolerance = 0.149842),"SpatialPolygonsDataFrame") 在这个answer 中完成，但是目前这实际上是如何工作的有点超出我的理解。但是，是的，我可以使用像素边界，事实上我已经这样做了，但这意味着创建 ID 列并合并两个数据帧，并且需要数百万个点的时间。我会为此发布另一个问题。

@stm4tt 使用您指向的答案我认为在这里不起作用，因为您的点网格未对齐。这个答案的关键是网格中心确实在 wgs lat long 但原始网格投影在另一个 crs 中。重新投影单元格进入原始 crs 使点对齐并适合SpatialPixels 转换。是否可以共享原始 NetCDF 数据以检查 crs ？

@Gilles 我明白了，所以我想唯一的方法是利用提供的像素角点，用它们构建多边形，构造一个sf 空间数据框并从那里开始（例如ggplot + geom_sf)。我会试一试。至于原始的 NetCDF，它是一个 600+MB 的文件，太大而无法共享（也不允许共享）。但我print(nc)ed 并粘贴了here。

我不确定你为什么想要多边形，但似乎还有很多工作流程可以将 R 中的 NetCDF 文件作为栅格（或类似栅格）直接读取。查看示例here 或here。

【参考方案1】：

这是执行此操作的一种方法。可能有一些更简单的方法，但这是可行的。

首先，我将使用 raster 包来操作坐标。我在这里创建的栅格是“非常规”的，因为它们包含的值是位置数据。但是使用栅格而不是矩阵可以访问一些有用的函数，例如extend，最有用的是resample，它的双线性插值函数将用于查找顶点

library(raster)
latr = raster(lat)
lonr = raster(lon)

find.vertices = function(m)
  r = raster(m)
  vertices = raster(matrix(nrow = dim(r)[1]+1, ncol = dim(r)[2]+1))
  x = extend(r, c(1,1))
  x[1,] = 2*x[2,] - x[3,]
  x[dim(x)[1],] = 2*x[dim(x)[1]-1,] - x[dim(x)[1]-2,]
  x[,1] = 2*x[,2] - x[,3]
  x[,dim(x)[2]] = 2*x[,dim(x)[2]-1] - x[,dim(x)[2]-2,]
  extent(vertices) = extent(r) + res(r)
  vertices = resample(x, vertices)


latv = find.vertices(lat)
lonv = find.vertices(lon)
df2 = data.frame(xc = lonv[], yc = latv[])

让我们绘制这些顶点以检查我们是否在轨道上：

ggplot(df, aes(y=lat, x=lon, fill=temp)) + 
  geom_tile(width=2, height=2) +
  geom_point(size=.1) +
  geom_point(aes(xc, yc), data=df2, inherit.aes =F) +
  borders('world', colour='gray50', size=.2) + 
  coord_quickmap(xlim=range(lon), ylim=range(lat)) +
  scale_fill_distiller(palette='Spectral') +
  theme_minimal()

现在我们从这些顶点创建一些Polygon

nx = NCOL(latv)
ny = NROW(lonv)
polys = list()
for (i in 1:length(data)) 
  x = col(data)[i]
  y = row(data)[i]
  polys[[i]] = Polygon(cbind(
      lonv[c((x-1)*ny + y, (x-1)*ny + y + 1, x*ny + y + 1, x*ny + y, (x-1)*ny + y)], 
      latv[c((x-1)*ny + y, (x-1)*ny + y + 1, x*ny + y + 1, x*ny + y, (x-1)*ny + y)]
    ))

将Polygon的列表转换为SpatialPolygonsDataFrame

Polys = sapply(1:length(polys), function(i) Polygons(polys[i], i))
SPolys = sapply(1:length(polys), function(i) SpatialPolygons(Polys[i], i))
SPolys = do.call(bind, SPolys)
SPolysdf = SpatialPolygonsDataFrame(SPolys, data.frame(data=as.vector(data)))

要在 ggplot 中绘制这个对象，我们需要转换为常规的data.frame。传统上，大多数人为此使用fortify。但是 ggplot 文档警告说这可能会被弃用，并建议使用 broom 包。我对扫帚不太熟悉，但我决定像这样遵循这个建议：

library(broom)
ggSPolysdf = tidy(SPolysdf)
ggSPolysdf = cbind(ggSPolysdf, data = rep(as.vector(data), each=5))

最后我们可以绘制：

ggplot(df, aes(y=lat, x=lon, fill=temp)) + 
  geom_polygon(aes(long,lat,fill=data, group = id), data=ggSPolysdf) +
  borders('world', colour='gray50', size=.2) + 
  coord_quickmap(xlim=range(lon), ylim=range(lat)) +
  scale_fill_distiller(palette='Spectral') +
  theme_minimal()

【讨论】：

优秀的答案。一个可能的改进是使用sf 而不是sp，我认为在未来的ggplot2 版本中也会支持这一点

【参考方案2】：

下面的解决方案基本上使用了来自@dww 的答案，并进行了一些似乎是获取数据所必需的更改（至少在我的平台上）。这些变化首先与定义上图中“倾斜像素”的多边形的定义有关；其次，关于如何将多边形压缩成数据框的问题。对于第二个问题，使用@SymbolixAU 建议的sf 包。

library(raster)
latr = raster(lat)
lonr = raster(lon)

find.vertices = function(m)
 r = raster(m)
 vertices = raster(matrix(nrow = dim(r)[1]+1, ncol = dim(r)[2]+1))
 x = extend(r, c(1,1))
 x[1,] = 2*x[2,] - x[3,]
 x[dim(x)[1],] = 2*x[dim(x)[1]-1,] - x[dim(x)[1]-2,]
 x[,1] = 2*x[,2] - x[,3]
 x[,dim(x)[2]] = 2*x[,dim(x)[2]-1] - x[,dim(x)[2]-2,]
 extent(vertices) = extent(r) + res(r)
 vertices = resample(x, vertices)


latv = find.vertices(lat)
lonv = find.vertices(lon)
df2 = data.frame(xc = lonv[], yc = latv[])

让我们绘制这些顶点以检查我们是否在轨道上：

ggplot(df, aes(y=lat, x=lon, fill=temp)) + 
  geom_tile(width=2, height=2) +
  geom_point(size=.1) +
  geom_point(aes(xc, yc), data=df2, inherit.aes=F) +
  borders('world', colour='gray50', size=.2) + 
  coord_quickmap(xlim=range(lon), ylim=range(lat)) +
  scale_fill_distiller(palette='Spectral') +
  theme_minimal()

现在我们从这些顶点创建多边形：

nx = NCOL(latv)

polys = list()
for (i in 1:length(data)) 
  x = col(data)[i]
  y = row(data)[i]

  polys[[i]] = Polygon(cbind(
      lonv[c((y-1)*nx + x, (y-1)*nx + x + 1, y*nx + x + 1, y*nx + x, (y-1)*nx + x)], 
      latv[c((y-1)*nx + x, (y-1)*nx + x + 1, y*nx + x + 1, y*nx + x, (y-1)*nx + x)]
    ))

将多边形列表转换为 SpatialPolygonsDataFrame：

Polys = sapply(1:length(polys), function(i) Polygons(polys[i], i))
SPolys = sapply(1:length(polys), function(i) SpatialPolygons(Polys[i], i))
SPolys = do.call(bind, SPolys)
SPolysdf = SpatialPolygonsDataFrame(SPolys, data.frame(data=as.vector(data)))

使用fortify 转换为数据框将通过以下两行来完成：（警告： 如@dww 所述，ggplot2 文档不建议使用此解决方案。 )

ggSPolysdf_0 = fortify(SPolysdf)
ggSPolysdf = cbind(ggSPolysdf_0, data = rep(as.vector(data), each=5))

另一种方法是使用sf 包。在下面的代码中，命令st_coordinates在ggplot2中扮演fortify的角色。请注意，使用本方法，变量名称在转换中丢失，需要手动恢复：

library(sf)
sfSPolys = st_as_sf(SPolysdf)
coord_xy_SPolys = st_coordinates(sfSPolys)
coord_xyz_SPolys = cbind(coord_xy_SPolys, data = rep(as.vector(data), each=5))
ggSPolysdf = as.data.frame(coord_xyz_SPolys)
colnames(ggSPolysdf) <- c("long", "lat", "piece", "id", "data")

最后我们可以绘制：

ggplot(df, aes(y=lat, x=lon, fill=temp)) + 
  geom_polygon(mapping=aes(long,lat,fill=data, group=id), data=ggSPolysdf) +
  borders('world', colour='gray50', size=.2) + 
  coord_quickmap(xlim=range(lon), ylim=range(lat)) +
  scale_fill_distiller(palette='Spectral') +
  theme_minimal()