案例分析：回归-克里金方法生成气温表面图（1）

,欢迎也在微信公众号查看。

背景知识

本部分内容主要是对气象站点的气温进行回归-克里金（regression-kriging，RK）插值生成气温表面图。RK方法是地统计中常用的方法，由趋势项和残差项构成。趋势项在该例中通过以气温为因变量，高程、站点到海岸线的距离为自变量，建立二着的线性模型，对残差建立克里金模型。使用RK的唯一要求是存在一个或多个与因变量显著相关的协变量。本例中将会用到七种回归算法，并利用10-fold CV来进行验证和均方根误差作为评估（Root-mean-square error，RMSE）。

10-fold CV(10-foldcross-validation),一种交叉验证方法，将数据集分成10份，轮流将其中9份用于训练，1份用于测试，循环10次，求准确度的平均值。

七种方法

Interpolation:
- 普通克里金(Ordinary Kriging,OK)
  - 相关参数说明（半变异函数、变程等）
  - 可查看以下视频https://www.bilibili.com/video/BV13E411T74X/
Regression:
- 广义线性模型(Generalized Linear Model,GLM)
- 广义加性模型 (Generalized Additive Models, GAM)
- 随机森林模型(Random Forest ,RF)
Regression-kriging:
- GLM + OK of residuals
- GAM + OK of residuals
- RF + OK of residuals

案例实现

本例中的数据是1950-2000年95个站点每天的平均气温，辅助数据为Elevation (Elev in meters a.s.l.), Distance to the coastline (distCoast in degrees); Latitude (Lat in degrees), and, Longitude (Lon in degrees).

数据下载

https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/CLIM_DATA_PT.zip

数据解压至F:/test/data-raw

library(raster)

## Loading required package: sp

setwd("F:/test/data-raw")
fl <- list.files("climData/rst", pattern = ".tif$", full.names = TRUE)
# 创建栅格数据集合
rst <- stack(fl)
# 修改名字
names(rst) <- c("distCoast", "Elev", "Lat", "Lon")
plot(rst)

读入每个站点的气温及对应的其它变量值

climDataPT<-read.csv("climData/clim_data_pt.csv")
climDataPT<-na.omit(climDataPT)
knitr::kable(head(climDataPT,n=10))

StationName	StationID	Lat	Lon	Elev	AvgTemp	distCoast
Sagres	1	36.98	-8.95	40	16.3	0.0000000
Faro	2	37.02	-7.97	8	17.0	0.0201246
Quarteira	3	37.07	-8.10	4	16.6	0.0090000
Vila do Bispo	4	37.08	-8.88	115	16.1	0.0360000
Praia da Rocha	5	37.12	-8.53	19	16.7	0.0000000
Tavira	6	37.12	-7.65	25	16.9	0.0458912
S. Br醩 de Alportel	7	37.17	-7.90	240	15.9	0.1853213
Vila Real Sto. Ant髇io	8	37.18	-7.42	7	17.1	0.0127279
Monchique	9	37.32	-8.55	465	15.0	0.1980000
Zambujeira	10	37.50	-8.75	106	15.0	0.0450000

将csv数据转为sp对象

proj4Str <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
statPoints <- SpatialPointsDataFrame(coords =climDataPT[,c("Lon","Lat")],
data= climDataPT,proj4string = CRS(proj4Str))
par(mfrow=c(1,2),mar=c(5,6,3,2))
plot(rst[["Elev"]], main="Elevation and weather stations",
xlab = "Longitude", ylab="Latitude")
plot(statPoints, add=TRUE)
hist(climDataPT$AvgTemp, xlab= "Temperature (ºC)", main="Annual avg. temperature")

从上面的左图可看出站点更多的分布在低纬沿海岸地区，右图反映出气温直方图是左偏的。接下来，我们通过计算相关系数矩阵来分析下气温与其它变量相关性强度。

library(corrplot)

corMat<-cor(climDataPT[,3:ncol(climDataPT)])
corrplot.mixed(corMat,number.cex=0.8,tl.cex=0.9,tl.col="black",
outline=F,mar = c(0,0,2,2),upper = "square",bg=NA)

相关系数图表明高程和纬度的相关性较大，与distCoast的相关性较小，在建立模型的过和中，是可考虑将其删除的。