案例分析:回归-克里金方法生成气温表面图(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的相关性 较小,在建立模型的过和中,是可考虑将其删除的。

已标记关键词 清除标记
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页