R_空间插值_必知必会（二）

欢迎关注天善智能，我们是专注于商业智能BI，人工智能AI，大数据分析与挖掘领域的垂直社区，学习，问答、求职一站式搞定！

对商业智能BI、大数据分析挖掘、机器学习，python，R等数据领域感兴趣的同学加微信：tstoutiao，邀请你进入数据爱好者交流群，数据爱好者们都在这儿。

作者：李誉辉  

四川大学在读研究生

前言

本文是R空间插值—必知必会的最后一篇，上一篇请戳：

R_空间插值_必知必会（一）

6、ggplot2绘图

6.1

rasterLayer转化为dataframe

 1library(raster)
 2library(sp)
 3library(dplyr)
 4library(magrittr)
 5
 6# 定义一个函数，将rasterLayer栅格数据转化为data.frame
 7#　将rasterLayer栅格数据转化为dataframe数据。
 8rasterL_to_DF <- function(climate_mask) {
 9  climate_mask_df <-  as.data.frame(cbind(coordinates(climate_mask), # 合并坐标数据和栅格值
10                                   values(climate_mask))
11                                   )  
12  climate_mask_df %<>% rename(climate_variable = V3) %<>% na.omit()
13  return(climate_mask_df)
14}
15
16# 调用函数
17TEM_mask_df <- rasterL_to_DF(climate_mask = TEM_mask) 
18
19# 生成等值线数据
20
21breaks_lines <- seq(min(TEM_1th$aver_TEM), max(TEM_1th$aver_TEM), length.out = 10)

6.2

绘图

 1library(ggplot2)
 2
 3ggplot_TEM <- function(temperature_data = TEM_mask_df) {
 4  ggplot(data = temperature_data) + 
 5    geom_raster(aes(x = x, y = y, fill = climate_variable)) + 
 6    # 增加等值线
 7    geom_contour(aes(x = x,y = y,z = climate_variable),
 8               color ="white", breaks=breaks_lines) +
 9    # 中国省级地图轮廓
10    geom_polygon(data = Chinaprovinces_df,
11                aes(x = long, y = lat, group = group), 
12                color = "black", fill = "transparent", size = 0.5) +
13    # 台湾地图
14    geom_polygon(data = Taiwan_df,
15                aes(x = long, y = lat, group = group), 
16                color = NA, fill = "grey", size = 0.5) + 
17    # 增加南海九段线
18    geom_line(data = Nine_lines, 
19              aes(x=long, y=lat, group=ID), 
20              colour="black", size=1) + 
21    # 各省名字，坐标为省会所在位置
22    geom_text(data = provinces,
23              aes(x = x, y = y, label = shortname), 
24              color = "green", size = 2) + 
25    coord_cartesian() + # geom_raster只能搭配笛卡尔坐标系
26    # 栅格颜色填充标度
27    scale_fill_gradient2(low = "blue", mid = "white", midpoint = 0, 
28                          high= "red",na.value = "grey",  # 设定缺失值为灰色
29                          name = "气温(℃)") +  #
30    labs(title = "中国平均气温分布图\n(2017年1月1日)", 
31        caption = "注：图中数据不含台湾地区") + 
32    theme_void()+
33    theme(
34      legend.position=c(0.2,0.2),
35      legend.background=element_blank(),
36      plot.title = element_text(colour = "magenta", size = 13, 
37                                face = "bold", hjust = 0.5),
38      legend.title = element_text(face = "bold", colour = "deeppink")
39    )
40}
41
42ggplot_TEM()
43

7、其他插值

7.1

多项式拟合

7.1.1 创建一个栅格产生函数

将：多项式模型、业务数据、空栅格、边界条件代入函数，即可产生栅格。

 1library(gstat)
 2library(sp)
 3library(raster)
 4
 5# 首先定义回归模型
 6
 7# 将自定义的多项式公式代入回归运算，然后将回归运算预测空栅格中的值，相当于插值计算。 
 8## 生成一个函数，可以直接调用多项式模型公式，和sp数据及空栅格数据
 9climate_mask_lm <- function(climate_sp, grd_climate, polynomial_function, boundary_sp) {
10
11  ### 添加变量X和Y
12  climate_sp$X <- coordinates(climate_sp)[,1]
13  climate_sp$Y <- coordinates(climate_sp)[,2]
14
15  ### 进行线性回归运算
16  lm_n <- lm(polynomial_function, data = climate_sp)
17
18  ### 将回归模型当做插值进行运算
19  climate_lm <- SpatialGridDataFrame(grd_climate, 
20                                data.frame(var1.pred = predict(lm_n, 
21                                                               newdata = grd_climate))) 
22
23  climate_raster   <- raster(climate_lm) # 栅格化
24  climate_mask <- mask(climate_raster, boundary_sp) # 筛选边界范围内的栅格
25
26  rm(lm_n, climate_lm, climate_raster)
27  return(climate_mask)
28}

7.1.2 一阶线性拟合

 1library(gstat)
 2library(sp)
 3library(raster)
 4library(tmap)
 5
 6# 自定义一个一阶线性公式
 7polynomial_1 <- as.formula(aver_TEM ~ X + Y) 
 8
 9# 调用上述自定义函数
10TEM_mask <- climate_mask_lm(climate_sp = TEM_sp, 
11                            grd_climate = grd_TEM,
12                            polynomial_function = polynomial_1,
13                            boundary_sp = Chinaboundary_noTaiwan_sp)
14
15# tmap绘图
16tmap_TEM()
17
18## ggplot2绘图
19TEM_mask_df <- rasterL_to_DF(climate_mask = TEM_mask) 
20ggplot_TEM()

7.1.3 二阶多项式拟合

 1library(gstat)
 2library(sp)
 3library(raster)
 4library(tmap)
 5
 6# 自定义一个二阶多项式公式
 7polynomial_2 <- as.formula(aver_TEM ~ X + Y + I(X^2)+I(Y^2) + I(X*Y))
 8
 9# 调用上述自定义函数
10TEM_mask <- climate_mask_lm(climate_sp = TEM_sp, 
11                            grd_climate = grd_TEM,
12                            polynomial_function = polynomial_2,
13                            boundary_sp = Chinaboundary_noTaiwan_sp)
14
15# tmap绘图
16tmap_TEM()
17
18## ggplot2绘图
19TEM_mask_df <- rasterL_to_DF(climate_mask = TEM_mask) 
20ggplot_TEM()

7.2

Kriging(克里金)插值

7.2.1 普通克里金插值

7.2.1.1 求拟合模型

求变异函数,首先绘制样本实验变异函数图(sample experimental variogram plot)。

1library(gstat)
2
3TEM_v <- variogram(aver_TEM ~ 1, data = TEM_sp, cloud = FALSE) # cloud = F只显示各个区间数字
4plot(TEM_v, plot.number = T) 
5

根据半变异图可知，已知点的自相关关系semivariance随着距离distance的增加而增加， 通过其分布，结合下图，可初步确定用线性函数或power函数来拟合，

拟合后绘图发现power函数更加恰当。

1library(gstat)
2
3TEM_v_fit <- fit.variogram(object = TEM_v,
4                           model = vgm(1, "Pow", 1))
5plot(TEM_v, TEM_v_fit) # 结果非常好

1TEM_v_fit <- fit.variogram(object = TEM_v,
2                           fit.ranges = FALSE, fit.sills = FALSE,
3                           model = vgm(psill = 18, model = "Sph", range = 28, nugget = 2.5))
4
5plot(TEM_v, TEM_v_fit) 
6

7.2.1.2 开始拟合运算

 1library(gstat)
 2library(raster)
 3
 4# 根据上面的拟合模型进行克里金插值的计算
 5TEM_krg <- gstat::krige(formula = aver_TEM ~ 1, 
 6                        model = TEM_v_fit,
 7                        locations = TEM_sp, # 数据点坐标
 8                        newdata = grd_TEM, # 需要插值点的位置
 9                        nmax = 15, nmin = 10 # 分布表示最多和最少搜索点的个数
10                        ) 
11
12TEM_raster <- raster(TEM_krg) # 栅格化
13TEM_mask <- mask(TEM_raster, Chinaboundary_noTaiwan_sp) # 筛选边界条件内的栅格
14
15rm(TEM_v, TEM_v_fit, TEM_krg, TEM_raster)

## [using ordinary kriging]

7.2.1.3 绘图

1library(tmap)
2library(ggplot2)
3
4# tmap绘图
5tmap_TEM()
6
7## ggplot2绘图
8TEM_mask_df <- rasterL_to_DF(climate_mask = TEM_mask) 
9ggplot_TEM()

7.2.2 泛克里金插值

用泛克里金法需谨慎，因其假定数据中存在覆盖趋势，
应该仅在了解数据中存在某种趋势并能够提供科学判断描述泛克里金法时，才可使用该方法
这在地质统计领域用得比较多，如矿藏的分布。
我这儿没有相关数据，仅仅用气温数据，可能不太准确，但是思路和流程是对的。
首先建立趋势模型，根据将大部分点，绘制样本实验变异函数图。
如下图，通过调节cutoff和width将大多数点显示在范围内，
可以使用plot.number = TRUE显示点的数量。

 1library(gstat)
 2
 3### 添加变量X和Y
 4TEM_sp2 <- TEM_sp # 复制一份
 5TEM_sp2$X <- coordinates(TEM_sp2)[,1]
 6TEM_sp2$Y <- coordinates(TEM_sp2)[,2]
 7trend_1 <- as.formula(aver_TEM ~ X + Y) 
 8
 9TEM_v <- variogram(object = trend_1, data = TEM_sp2, cloud = FALSE,
10                   cutoff = 30, # cutoff为对角线长度调整，其与width会相互影响
11                   width = 2) # width表示相邻2个点之间的distance，width越小，点越多
12plot(TEM_v) 

1TEM_v_fit <- fit.variogram(object = TEM_v,
2                           fit.ranges = FALSE, fit.sills = FALSE,
3                           model = vgm(psill = 18, model = "Sph", range = 28, nugget = 2.5))
4
5plot(TEM_v, TEM_v_fit) 
6

运算量非常大，经常溢出，这里不进行运算。

 1library(gstat)
 2library(raster)
 3
 4# 根据上面的拟合模型进行克里金插值的计算
 5TEM_krg <- gstat::krige(formula = trend_1,
 6                        locations = TEM_sp2, # 数据点坐标
 7                        newdata = grd_TEM, # 需要插值点的位置
 8                        model = TEM_v_fit
 9                        ) 
10
11TEM_raster <- raster(TEM_krg) # 栅格化
12TEM_mask <- mask(TEM_raster, Chinaboundary_noTaiwan_sp) # 筛选边界条件内的栅格
13
14rm(TEM_v, TEM_v_fit, TEM_krg, TEM_raster, TEM_sp2)

1library(tmap)
2library(ggplot2)
3
4# tmap绘图
5tmap_TEM()
6
7## ggplot2绘图
8TEM_mask_df <- rasterL_to_DF(climate_mask = TEM_mask) 
9ggplot_TEM()

7.3

akima插值

akima插值不支持sp数据对象的插值，只支持dataframe和matrix对象插值。

插值结果为dataframe对象, 只能形成SpatialPixelsDataFrame栅格类型，
与前面的sp对象插值不同，sp对象插值结果为SpatialGridDataFrame栅格类型。
目前dataframe对象插值更简单，
但是SpatialPixelsDataFrame对象不支持多个多边形边界进行筛选。 over()不支持多个多边形边界进行筛选 ，mask()不支持SpatialPixelsDataFrame对象。
所以只能一个个边界进行筛选，然后索引内部元素进行合并。

7.3.1 插值运算

 1library(akima)
 2library(sp)
 3library(raster)
 4
 5# 对整个TEM_1th进行插值
 6TEM_interp <- interp(x = TEM_1th$long, y = TEM_1th$lat, z = TEM_1th$aver_TEM,
 7                     xo = seq(60, 140, by = 0.1), # 指定插值经度范围
 8                     yo = seq(10, 60, 0.1),       # 指定插值纬度范围
 9                     linear = FALSE, # 表示是线性插值还是样条插值
10                     extrap = TRUE) # 表示是否接受外插，有的栅格只能外插才能得到，
11
12# 生成网格数据
13TEM_grd <- SpatialPoints(expand.grid(x=TEM_interp$x, y = TEM_interp$y))
14TEM_grd <- SpatialPixelsDataFrame(TEM_grd, data.frame(kde = array(TEM_interp$z,
15                                                            length(TEM_interp$z))))
16# 分离地图边界
17df_as_sp <- function(map_df, area) { # x,y指定经纬度
18  map_subset <-  subset(map_df, AREA == area)
19  Sr1 <- Polygon(cbind(map_subset$long, map_subset$lat)) 
20  Srs1 <- Polygons(list(Sr1), ID = "1") 
21  SpP <- SpatialPolygons(Srl = list(Srs1), 1:1)
22  partmap_sp <- SpatialPolygonsDataFrame(
23    Sr = SpP,
24    data = data.frame(Names = "coords", row.names = row.names(SpP)))
25  return(partmap_sp)
26}
27
28Mainland_sp <- df_as_sp(Chinaboundary_df, 954.943)  
29Hainan_sp <- df_as_sp(Chinaboundary_df, 2.903) 
30
31# 筛选各个边界内的栅格数据
32Mainland_overcheck <- !is.na(sp::over(x = TEM_grd, y = Mainland_sp)) 
33Hainan_overcheck <- !is.na(sp::over(x = TEM_grd, y = Hainan_sp))
34Mainland_grd <- TEM_grd[Mainland_overcheck[,1], ]
35Hainan_grd <- TEM_grd[Hainan_overcheck[,1], ]
36
37# 栅格合并
38Mainland_grd <-cbind(Mainland_grd@coords,Mainland_grd@data) # 
39Hainan_grd <- cbind(Hainan_grd@coords,Hainan_grd@data)
40
41grd_bind_noTaiwan <- rbind(Mainland_grd, Hainan_grd)
42
43rm(TEM_interp, TEM_grd, df_as_sp,
44   Mainland_overcheck, Hainan_overcheck, 
45   Mainland_grd, Hainan_grd )# 移除中途数据
46
47# 生成等值线数据
48breaks_lines <- seq(min(TEM_1th$aver_TEM), max(TEM_1th$aver_TEM), length.out = 10)

7.3.2 ggplot2绘图

 1library(ggplot2)
 2
 3ggplot(data = grd_bind_noTaiwan) + 
 4  # 所有栅格
 5  geom_raster(aes(x=x,y=y,fill=kde)) +    
 6  # 增加等值线
 7  geom_contour(aes(x=x,y=y,z=kde),
 8               color ="white",breaks=breaks_lines) + 
 9  # 中国省级地图轮廓
10  geom_polygon(data = Chinaprovinces_df,
11               aes(x = long, y = lat, group = group), 
12               color = "black", fill = "transparent", size = 0.5) + 
13  # 台湾地图
14  geom_polygon(data = Taiwan_df,
15               aes(x = long, y = lat, group = group), 
16               color = NA, fill = "grey", size = 0.5) + 
17  # 各省名字，坐标为省会所在位置
18  geom_text(data = provinces,
19            aes(x = x, y = y, label = shortname), 
20            color = "green", size = 2) + 
21  # 增加南海九段线
22  geom_line(data = Nine_lines, 
23            aes(x=long, y=lat, group=ID), 
24            colour="black", size=1) + 
25  coord_cartesian() + # geom_raster只能搭配笛卡尔坐标系
26   # 栅格颜色填充标度
27  scale_fill_gradient2(low = "blue", mid = "white", midpoint = 0, 
28                       high= "red",na.value = "grey",  # 设定缺失值为灰色
29                       name = "气温(℃)") +  #
30  labs(title = "中国平均气温分布图\n(2017年1月1日)", 
31       caption = "注：图中数据不含台湾地区") + 
32  theme_void()+
33  theme(
34    legend.position=c(0.2,0.2),
35    legend.background=element_blank(),
36    plot.title = element_text(colour = "magenta", size = 13, 
37                              face = "bold", hjust = 0.5),
38    legend.title = element_text(face = "bold", colour = "deeppink")
39  )

8、插入小地图

tmap中只能插入与leaflet中类似的小地图，即在线小地图。
实际上还是用PS截图拼图更加方便，甚至用PPT也行。
小地图会增加上百行代码。

参

考来源

R语言空间插值的几种方法及案例应用

gstat插值参数确定

https://mgimond.github.io/Spatial/spatial-interpolation.html

tmap空间插值

https://mgimond.github.io/Spatial/interpolation-in-r.html

R中的点插值

https://www.cdrc.ac.uk/wp-content/uploads/2016/11/Practical_11.html

tmap实例

https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html

tmap绘制人口调查地图

http://zevross.com/blog/2018/10/02/creating-beautiful-demographic-maps-in-r-with-the-tidycensus-and-tmap-packages/

tmap分面及布局

https://geocompr.robinlovelace.net/adv-map.html

sf介绍

https://cran.r-project.org/web/packages/sf/vignettes/sf1.html

SpatialLinesDataFrame创建

https://gis.stackexchange.com/questions/163286/how-do-i-create-a-spatiallinesdataframe-from-a-dataframe

SpatialPolygonsDataFrame创建

https://www.rdocumentation.org/packages/sp/versions/1.3-1/topics/SpatialPolygonsDataFrame-class

SPDF转化为SLDF

https://gis.stackexchange.com/questions/200679/convert-spatialpolygonsdf-boundaries-to-spatiallinesdf-keeping-information-on-po

sp对象及栅格数据介绍(最全)

https://rspatial.org/spatial/8-rastermanip.html

合并多个SPDF

https://gis.stackexchange.com/questions/155328/merging-multiple-spatialpolygondataframes-into-1-spdf-in-r

Kriging插值

https://rpubs.com/nabilabd/118172

IDW插值

https://nceas.github.io/oss-lessons/spatial-data-gis-law/4-tues-spatial-analysis-in-r.html

akima插值与ggplot2绘图

https://stackoverflow.com/questions/50533738/best-method-of-spatial-interpolation-for-geographic-heat-contour-maps

IDW插值与ggplot2

http://aasa.ut.ee/LOOM02331/R_idw_interpolation.html

kknn插值与ggplot2

https://timogrossenbacher.ch/2018/03/categorical-spatial-interpolation-with-r/

反距离加权插值

https://www.jianshu.com/p/b38c5e464d16

将rasterLayer对象转化为SpatialGrid/SpatialPixels对象

https://gis.stackexchange.com/questions/43707/how-to-produce-spatial-grid-from-raster

KNN插值原理

http://www.kuqin.com/algorithm/20120817/329048.html

往期精彩：

• R_空间插值_必知必

  会（一）

• ggplot2图集汇总（一）

• R_ggplot2地理信息可视化_史上最全（一）

• R语言中文社区2018年终文章整理（作者篇）

• R语言中文社区2018年终文章整理（类型篇）

公众号后台回复关键字即可学习

回复 爬虫            爬虫三大案例实战
回复 Python       1小时破冰入门
回复 数据挖掘     R语言入门及数据挖掘
回复 人工智能     三个月入门人工智能
回复 数据分析师  数据分析师成长之路 
回复 机器学习     机器学习的商业应用
回复 数据科学     数据科学实战
回复 常用算法     常用数据挖掘算法