在R中创建用于克里金插值等插值方法的数据对象

Question

我有不同位置（X1、X2等）的每日平均温度数据，并希望用它们插值出地图。我通过从格式化的Excel表格中加载它们创建了一个长格式数据对象，例如：

library(reshape2)
tempdata <- read.csv("...", sep=";")
names(tempdata) <- c("date", paste("X", 1:73))
head(tempdata)
#    date  X1  X2  X3  X4  X5  X6  X7
# 1    1  7.3 6.6 6.7 5.8 6.1 6.1 5.5
# 2    2  7.5 6.6 6.6 5.6 4.8 4.7 3.9
# 3    3  8.8 7.7 7.6 7.0 7.0 6.0 5.8
# 4    4  8.5 7.4 7.5 7.0 7.3 5.9 5.5
# 5    5  7.7 6.7 6.9 6.1 6.8 5.1 4.1
# 6    6  7.5 6.7 6.8 6.0 6.4 5.0 4.1

各位置 X1、X2 等的纬度和经度相同：

lat.lon <- read.csv("...", sep=";")
rownames(lat.lon) <- c(paste0("X",1:73))
head(lat.lon)
#     latitude longitude
#  X1  54.1650    6.3458
#  X2  54.1667    7.4500
#  X3  54.1832    7.8856
#  X4  55.0114    8.4158
#  X5  54.5068    9.5393
#  X6  54.5214   11.0522

我把它们合并成了长格式：

res <- merge(
  melt(tempdata, id.vars="date"), 
  lat.lon, 
  by.x="variable", by.y="row.names"
)
head(res)
#  variable     date value latitude longitude
#       X1        1   9.9   54.165    6.3458
#       X1        2   8.9   54.165    6.3458
#       X1        3   7.8   54.165    6.3458
#       X1        4   9.2   54.165    6.3458
#       X1        5   8.7   54.165    6.3458
#       X1        6   8.4   54.165    6.3458

使用

coordinates(res) = ~longitude+latitude

我可以使用spplot在正确的位置上绘制它们，同时还包括国家边界：

library(maptools)
load(url('http://gadm.org/data/rda/DEU_adm0.RData'))
GE <- gadm
GE <- spChFIDs(GE, paste("GE", rownames(GE), sep = "_"))
spplot(res["value"], sp.layout = list("sp.polygons", GE), col.regions=bpy.colors(20))

我想在观测的单个日期使用IDW方法，但我找到的包中的idw方法（例如gstat）需要其他“格网”数据对象。我该如何创建这样的数据对象以便使用这些方法进行插值？

- MojioMS

Answer 1

像这样

加载必要的包

（注：原文已经是中文，因此无需翻译）

kpacks <- c('sp','rgdal', 'gstat', 'raster')
new.packs <- kpacks[!(kpacks %in% installed.packages()[,"Package"])]
if(length(new.packs)) install.packages(new.packs)
lapply(kpacks, require, character.only=T)
remove(kpacks, new.packs)

数据（wrld_simpl）

一个用于处理的投影坐标系

p.utm33n <- CRS("+init=epsg:32633") # UTM 33N Landsat Images

一个国家（我特别喜欢这个国家）

ago <- wrld_simpl[wrld_simpl@data$NAME == 'Angola',]

将其投影到UTM 33S。

ago <- spTransform(ago, p.utm33n)

在多边形内随机采样一些点

ago_p <- spsample(ago, type="random", n=25)    
plot(ago, col = 'grey' , axes = T)
plot(ago_p, add = T)

三天的虚构温度数据

tdata <- data.frame(x=rep(coordinates(ago_p)[,1], 3), 
                    y=rep(coordinates(ago_p)[,2], 3),
                    temp=runif(75, 12,35),
                    day = rep(1:3, each = 25))

将其转换为spatialPointDataFrame对象进行管理

coordinates(tdata) <- ~x+y 

proj4string(tdata) <- CRS(proj4string(ago))

由于我不知道你的基础地图，所以我将使用我上面选择的国家。 基础图层必须是一个SpatialPixelDataBase。我将使用一个rasterLayer进行操作。

rago <- raster(extent(ago))
res(rago) <- c(10000,10000)
rago[] <- 1
proj4string(rago) <- CRS(proj4string(ago))
r_ago <- mask(rago, ago)
#plot(r_ago)
grid_ago <- as(r_ago, 'SpatialPointsDataFrame')
grid_ago <- grid_ago[!is.na(grid_ago@data$layer), ]
gridded(grid_ago) <- TRUE

我现在可以运行来自gstat的idw()。我将使用day == 1的数据运行。

idw_ago <- idw(temp ~ 1, tdata[tdata$day == 1, ], grid_ago, idp = 2.5)

最后将其绘制出来

spplot(idw_ago, "var1.pred")

现在有了你的数据，这是我在你的问题中缺失的。同样的方法。

library(latticeExtra)
p.dutch <- CRS("+init=epsg:28991") # Dutch National Grid EPSG:28991
load(url('http://gadm.org/data/rda/DEU_adm0.RData'))
ger <- gadm
ger <- spChFIDs(ger, paste("ger", rownames(ger), sep = "_"))
ger <- spTransform(ger, p.dutch)
ger_p <- spsample(ger, type="random", n=25)
plot(ger, col = 'yellow', border = NA, axes = T, cex.axis = 0.6)
plot(ger_p, add = T, pch = 20)

tdata <- data.frame(x=rep(coordinates(ger_p)[,1], 3), 
                    y=rep(coordinates(ger_p)[,2], 3),
                    temp=runif(75, 12,35),
                    day = rep(1:3, each = 25))    
coordinates(tdata) <- ~x+y 
proj4string(tdata) <- CRS(proj4string(ger))    
rger <- raster(extent(ger))
res(rger) <- c(10000,10000)
rger[] <- 1
proj4string(rger) <- CRS(proj4string(ger))
r_ger <- mask(rger, ger)
plot(r_ger)
grid_ger <- as(r_ger, 'SpatialPointsDataFrame')
grid_ger <- grid_ger[!is.na(grid_ger@data$layer), ]
gridded(grid_ger) <- TRUE
idw_ger <- idw(temp ~ 1, tdata[tdata$day == 1, ], grid_ger, idp = 2.5)
spplot(idw_ger, "var1.pred") +
latticeExtra::layer(sp.polygons(ger, fill = NA, col = 'blue')) +
latticeExtra::layer(sp.points(tdata[tdata$day == 1, ],
                                fill = NA, col = 'red'))

希望它能有所帮助

- Paulo E. Cardoso

---

@MojioMS 我需要更多的信息。您是否使用完全相同的数据或提供其他数据集？ - Paulo E. Cardoso

上述代码对我来说完全正常，不需要任何更改。 - Paulo E. Cardoso

我终于让它工作了，但我不得不稍微改变我的数据框架。感谢您的帮助。 - MojioMS