我有不同位置(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)需要其他"网格化"数据对象.我如何创建这样的数据对象,以便用这些方法插入它们?
像这样的东西
加载必要的包
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)
一些假想的温度数据为3天
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
我现在可以idw()
从gstat 运行了.我将使用日期== 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'))
希望能帮助到你