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

作者： 李誉辉

四川大学在读研究生

前言

本文是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年终文章整理（类型篇）

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

本文由R语言中文社区 创作，采用 知识共享署名-相同方式共享 3.0 中国大陆许可协议 进行许可。

转载、引用前需联系作者，并署名作者且注明文章出处。

本站文章版权归原作者及原出处所有 。内容为作者个人观点， 并不代表本站赞同其观点和对其真实性负责。本站是一个个人学习交流的平台，并不用于任何商业目的，如果有任何问题，请及时联系我们，我们将根据著作权人的要求，立即更正或者删除有关内容。本站拥有对此声明的最终解释权。