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

Question

我已经将此问题作为在不规则网格问题上绘制数据的有效方法的一部分,但一般反馈是将原始问题拆分为更易于管理的块.因此,这个新问题.

我使用组织在不规则二维网格上的卫星数据,其尺寸为扫描线(沿轨道尺寸,即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包裹？

我可以看到它在这个问题的答案中完成,但是sf对于我来说,使用的答案有点不清楚,因为它涉及不同的投影和可能的常规网格,而在我的情况下,我想我不必重新投射我的数据,而且,我的数据不在常规网格上.

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

Answer 1

这是一种方法.可能有一些更简单的东西,但这是有效的.

首先,我将使用栅格包来操纵坐标.我在这里创建的栅格是"非常规的",因为它们包含的值是位置数据.但是使用栅格而不是矩阵可以访问一些有用的函数,例如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文档警告说这可能会被弃用,建议改用扫帚包.我对扫帚不太熟悉,但我决定这样遵循这样的建议:

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