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

欢迎关注天山情报。我们是一个专注于商业智能BI、人工智能AI和大数据分析与挖掘的垂直社区。学习、问答、求职一站搞定！对商业智能BI、大数据分析与挖掘、机器学习、python、R等数据领域感兴趣的同学加微信:tstoutiao，邀请你加入数据爱好者交流群，数据爱好者都在这里。

四川大学是研究生

前言

本文是最后一篇关于R空之间插值的文章，请戳前一篇:

R_空之间的插值_必须知道且必须知道(I)

6.ggplot2图纸

6.1

光栅层被转换成数据帧

1图书馆(光栅)

2图书馆(sp)

3库(dplyr)

4 library(magritter)

五

6#定义了将栅格图层栅格数据转换为数据的函数

7#将栅格图层栅格数据转换为数据帧数据。

8rasterL _ to _ DF & lt-函数(climate_mask) {

9 climate _ mask _ df & lt-as.data.frame (cbind(坐标(climate _ mask)，#组合坐标数据和栅格值

10个值(气候掩码))

11 )

12 climate _ mask _ df % & lt>。% rename(climate _ variable = V3)% & lt；>。% na .省略()

13返回(climate_mask_df)

14}

15

16#调用函数

17TEM _ mask _ df & lt-rasterL _ to _ DF(climate _ mask = TEM _ mask)

18

19#生成等值线数据

20

21breaks _ lines & lt。- seq(min(TEM_1th$aver_TEM)，max(TEM_1th$aver_TEM)，length.out = 10)

6.2

草稿

1图书馆(gplot2)

2

3ggplot _ TEM & lt-函数(温度_数据= TEM_mask_df) {

4g plot(数据=温度_数据)+

5 geom_raster(aes(x = x，y = y，fill = climate_variable)) +

将等值线添加到6 #

7 geom _等高线(aes(x = x，y = y，z = climate_variable)、

8 color ="white "，breaks = breaks _ lines

9 #中国省级地图概要

10 geom _多边形(数据=中国省份_df，

11个aes(x = long，y = lat，group = group)，

12 color = "黑色"，fill = "透明"，size = 0.5) +

13 #台湾地图

14 geom_polygon(数据= Taiwan_df，

15个aes(x =长，y = lat，group = group)，

16色= NA，填充=“灰色”，尺寸= 0.5) +

17 #增加南海九段线

18 geom_line(数据=九条线，

19个AES(x =长，y=lat，group=ID)，

20 colour = " black "，size=1) +

21 #省名，坐标是省会

22 geom_text(数据=省份，

23个aes(x = x，y = y，label = shortname)，

24 color = "green "，size = 2) +

25 coord _ cartesian()+# geom _ raster只能与笛卡尔坐标系匹配

26 #网格颜色填充比例

27 scale_fill_gradient2(low = "蓝色"，mid = "白色"，中点= 0，

28 high= "红色"，na.value = "灰色"，#将缺少的值设置为灰色

29 name = "温度(℃)"+#

30个实验室(标题为“中国平均温度分布图 n(2017年1月1日)”，

31题注=“注:图中数据不包括台湾省地区”)+

32主题_void()+

33主题(

34图例.位置=c(0.2，0.2)，

35 legend . background = element _ blank()，

36 plot . title = element _ text(color = "品红"，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创建网格生成功能

通过将多项式模型、业务数据、空网格和边界条件代入函数，可以生成网格。

1图书馆(gstat)

2图书馆(sp)

3库(栅格)

四

5#首先定义回归模型

六

7#将自定义多项式公式代入回归运算，再通过回归运算预测空网格中的值，相当于插值计算。

8##生成一个可以直接调用多项式模型公式、sp数据和空栅格数据的函数

9climate _ mask _ lm & lt-函数(climate_sp，grd_climate，多项式_function，boundary_sp) {

10

11 ###添加变量x和y

12 climate _ sp $ X & lt-坐标(climate_sp)[，1]

13气候_标准普尔$ Y & lt-坐标(climate_sp)[，2]

14

15 ###执行线性回归运算

16 lm _ n & lt- lm(多项式函数，数据= climate_sp)

17

18 ###将回归模型作为插值操作

19气候_ lm & lt-空间网格框架(grd_climate，

20 data . frame(var 1 . pred = predict(lm _ n，

21 newdata = grd_climate)))

22

23气候_栅格<。-光栅(climate_lm) #光栅化

24气候面具<。-遮罩(climate _ raster，boundary _ sp) #过滤边界内的栅格

25

26 rm(lm_n，climate_lm，climate_raster)

27返回(气候面具)

28}

7.1.2一阶线性拟合

1图书馆(gstat)

2图书馆(sp)

3库(栅格)

4图书馆(tmap)

五

6#定制一阶线性公式

7多项式_ 1 & lt- as.formula(aver_TEM ~ X + Y)

八

9#调用上面的自定义函数

10TEM _ mask & lt-climate _ mask _ lm(climate _ sp = TEM _ sp，

11 grd_climate = grd_TEM，

12多项式_函数=多项式_1，

13 boundary _ sp = China boundary _ notiwan _ sp)

14

15# tmap图纸

16tmap _ TEM()

17

18## ggplot2图纸

19TEM _ mask _ df & lt-rasterL _ to _ DF(climate _ mask = TEM _ mask)

20ggplot_TEM()

7.1.3二阶多项式拟合

1图书馆(gstat)

2图书馆(sp)

3库(栅格)

4图书馆(tmap)

五

6#自定义二阶多项式公式

7多项式_ 2 & lt-as . formula(aver _ tem ~ x+y+i(x^2)+i(y^2)+I(x * y))

八

9#调用上面的自定义函数

10TEM _ mask & lt-climate _ mask _ lm(climate _ sp = TEM _ sp，

11 grd_climate = grd_TEM，

12多项式_函数=多项式_2，

13 boundary _ sp = China boundary _ notiwan _ sp)

14

15# tmap图纸

16tmap _ TEM()

17

18## ggplot2图纸

19TEM _ mask _ df & lt-rasterL _ to _ DF(climate _ mask = TEM _ mask)

20ggplot_TEM()

7.2

克里金插值

7.2.1普通克里金插值

7.2.1.1拟合模型

要找到变差函数，首先绘制样本实验变差图。

1图书馆(gstat)

2

3TEM _ v & lt-变异函数(aver _ tem ~ 1，data = tem _ sp，cloud = false) # cloud = f仅显示每个间隔的数字

4批次(TEM_v，绘图.编号= T)

五

根据半方差图，已知点的自相关关系随着距离的增加而增加。通过其分布，结合下图，可以初步确定用线性函数或幂函数拟合。

拟合后发现幂函数更合适。

1图书馆(gstat)

2

3TEM _ v _ fit & lt-拟合变差函数(object = TEM_v，

4型号= vgm(1，“Pow”，1))

5plot(TEM_v，TEM_v_fit) #的结果很好

1TEM _ v _ fit & lt-拟合变差函数(object = TEM_v，

2拟合范围=假，拟合门槛=假，

3型号= vgm(psill = 18，model = "Sph "，范围= 28，熔核= 2.5))

四

5批次(TEM_v，TEM_v_fit)

六

7.2.1.2开始试衣手术

1图书馆(gstat)

2图书馆(光栅)

三

4#根据上述拟合模型计算克里金插值

5TEM _ krg & lt- gstat::krige(公式= aver_TEM ~ 1，

6型号= TEM_v_fit，

7个位置= TEM_sp，#数据点坐标

8 newdata = grd_TEM，#需要插值点的位置

9 nmax = 15，nmin = 10 #分布表示最大和最小搜索点数

10 )

11

12TEM _ raster & lt-光栅(TEM_krg) #光栅化

13TEM _ mask & lt-遮罩(tem _ raster，China boundary _ notaiwan _ sp) #在边界条件内过滤栅格

14

15rm(TEM_v，TEM_v_fit，TEM_krg，TEM_raster)

# #[使用普通克里金法]

7.2.1.3绘画

1图书馆(tmap)

2图书馆(gplot2)

三

4# tmap图纸

5tmap _ TEM()

六

7## ggplot2图纸

8TEM _ mask _ df & lt-rasterL _ to _ DF(climate _ mask = TEM _ mask)

9ggplot_TEM()

7.2.2泛克里金插值

泛克里金法应该谨慎，因为它假设数据中存在覆盖趋势。

只有当我们知道数据中存在一定的趋势，并且能够提供科学的判断来描述泛克里金法时，才应该使用这种方法

这在地质统计学领域有着广泛的应用，如矿床分布。

我这里没有相关数据。只使用温度数据可能不太准确，但思路和过程是正确的。

首先建立趋势模型，根据大部分点绘制样本实验的变差图。

如下所示，通过调整截止值和宽度，大多数点都显示在该范围内。

可以使用plot.number = TRUE来显示点数。

1图书馆(gstat)

2

3###添加变量x和y

4TEM _ sp2 & lt-TEM_sp #制作副本

5TEM _ sp2 $ X & lt-坐标(TEM_sp2)[，1]

6TEM _ sp2 $ Y & lt-坐标(TEM_sp2)[，2]

7trend_1 <。- as.formula(aver_TEM ~ X + Y)

八

9TEM _ v & lt-变差函数(object = trend_1，data = TEM_sp2，cloud = FALSE，

10截断= 30，#截断是对角线长度的调整，这将与宽度相互影响

11 width = 2) # width表示两个相邻点之间的距离，宽度越小，点越多

12批次(TEM_v)

然后根据点的趋势和对照表，选择合适的拟合模型。不同的型号在vgm()中有不同的参数。

这里我们选择球形模型。其三个参数分别代表以下含义:

1TEM _ v _ fit & lt-拟合变差函数(object = TEM_v，

2拟合范围=假，拟合门槛=假，

3型号= vgm(psill = 18，model = "Sph "，范围= 28，熔核= 2.5))

四

5批次(TEM_v，TEM_v_fit)

六

操作量很大，经常溢出，这里不做操作。

1图书馆(gstat)

2图书馆(光栅)

三

4#根据上述拟合模型计算克里金插值

5TEM _ krg & lt- gstat::krige(formula = trend_1，

6个位置= TEM_sp2，#数据点坐标

7 newdata = grd_TEM，#需要插值点的位置

8型号= TEM_v_fit

9 )

10

11TEM _ raster & lt-光栅(TEM_krg) #光栅化

12TEM _ mask & lt-遮罩(tem _ raster，China boundary _ notaiwan _ sp) #在边界条件内过滤栅格

13

14毫米(TEM_v，TEM_v_fit，TEM_krg，TEM_raster，TEM_sp2)

1图书馆(tmap)

2图书馆(gplot2)

三

4# tmap图纸

5tmap _ TEM()

六

7## ggplot2图纸

8TEM _ mask _ df & lt-rasterL _ to _ DF(climate _ mask = TEM _ mask)

9ggplot_TEM()

7.3

阿基马插值

Akima插值不支持sp数据对象的插值，只支持数据帧和矩阵对象的插值。

插值的结果是一个数据帧对象，只能形成SpatialPixelsDataFrame栅格类型。

与之前的sp对象插值不同，sp对象插值结果是SpatialGridDataFrame栅格类型。

目前，数据帧对象的插值比较简单，

但是，SpatialPixelsDataFrame对象不支持多个多边形边界进行过滤。Over()不支持多个多边形边界进行过滤，mask()不支持SpatialPixelsDataFrame对象。

因此，只能过滤一个边界，然后可以索引内部元素进行合并。

7.3.1插值操作

1图书馆(akima)

2图书馆(sp)

3库(栅格)

四

# 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线性=假，#表示是线性插值还是样条插值

10 extrap = TRUE) #表示是否接受外推，有些网格只能通过外推得到。

11

12#生成网格数据

13TEM _ grd & lt-Spatial Points(expand . grid(x = TEM _ interp $ x，y = TEM_interp$y))

14TEM_grd <。-SpatialPixelsDataFrame(TEM _ grd，data . frame(kde = array(TEM _ interp $ z，

15长度(TEM_interp$z)))

16#分隔地图边界

17df_as_sp <。-函数(map_df，area) {# x，y指定纬度和经度

18 map _ subset & lt-子集(map_df，AREA == area)

19 Sr1 <。-多边形(cbind(map_subset$long，map_subset$lat))

20 Srs1 <。-多边形(列表(Sr1)，ID = "1 ")

21 SpP <。-空间多边形(Srl =列表(Srs1)，1:1)

22 partmap_sp <。- SpatialPolygonsDataFrame(

23 Sr = SpP，

24 data = data . frame(Names = " coords "，row.names = row.names(SpP)))

25返回(partmap_sp)

26}

27

28Mainland _ sp & lt- df_as_sp(Chinaboundary_df，954.943)

29海南_ sp & lt- df_as_sp(Chinaboundary_df，2.903)

30

31#过滤每个边界内的栅格数据

32Mainland _ overcheck & lt- !is.na(sp::over(x = TEM_grd，y =大陆_sp))

33海南_过度检查& lt- !is.na(sp::over(x = TEM_grd，y =海南_sp))

34Mainland _ grd & lt。-TEM _ grd[大陆_过检查[，1]，]

35海南_ grd & lt-TEM _ grd[海南_overcheck[，1]，]

36

37#栅格合并

38Mainland _ grd & lt。-cbind(大陆_grd@coords，大陆_grd@data) #

39海南_ grd & lt- cbind(海南_grd@coords，海南_grd@data)

40

41grd _ bind _ noTaiwan & lt-rbnd(大陆_grd，海南_grd)

四十二岁

43rm(TEM_interp，TEM_grd，df_as_sp，

44大陆_超检，海南_超检，

45大陆_grd，海南_grd )#删除中途数据

46

47#生成等值线数据

48breaks _ lines & lt- seq(min(TEM_1th$aver_TEM)，max(TEM_1th$aver_TEM)，length.out = 10)

7 . 3 . 2 ggplot 2图纸

1图书馆(gplot2)

2

3gg plot(data = grd _ bind _ NotiWan)+

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 _多边形(数据=中国省份_df，

11个aes(x = long，y = lat，group = group)，

12 color = "黑色"，fill = "透明"，size = 0.5) +

13 #台湾地图

14 geom_polygon(数据= Taiwan_df，

15个aes(x =长，y = lat，group = group)，

16色= NA，填充=“灰色”，尺寸= 0.5) +

17 #省名，坐标是省会的位置

18 geom_text(数据=省份，

19个aes(x = x，y = y，label = shortname)，

20 color = "green "，size = 2) +

21 #增加南海九段线

22 geom_line(数据=九条线，

23个AES(x =长，y=lat，group=ID)，

24 colour = " black "，size=1) +

25 coord _ cartesian()+# geom _ raster只能与笛卡尔坐标系匹配

26 #网格颜色填充比例

27 scale_fill_gradient2(low = "蓝色"，mid = "白色"，中点= 0，

28 high= "红色"，na.value = "灰色"，#将缺少的值设置为灰色

29 name = "温度(℃)"+#

30个实验室(标题为“中国平均温度分布图 n(2017年1月1日)”，

31题注=“注:图中数据不包括台湾省地区”)+

32主题_void()+

33主题(

34图例.位置=c(0.2，0.2)，

35 legend . background = element _ blank()，

36 plot . title = element _ text(color = "品红"，size = 13，

37 face = "bold "，hjust = 0.5)，

38 legend . title = element _ text(face = " bold "，colour = " deeppink ")

39 )

8.插入一张小地图

只有类似于传单中的迷你地图才能插入tmap，即在线迷你地图。

其实用PS截图拼图更方便，甚至是PPT。

小地图加几百行代码。

参加

测试源

r语言空之间插值的几种方法及实例应用

gstat插值参数的确定

https://mgi mond . 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/Practual _ 11 . html

tmap示例

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

Tmap绘制普查地图

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

Tmap面和布局

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

sf简介

https://cran . r-project . org/web/packages/SF/晕映/sf1.html

空间线数据帧创建

https://GIS . stackexchange . com/questions/163286/how-do-I-create-a-spatialinedataframe-from-a-data frame

空间多边形数据帧创建

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

SPDF变成了SLDF

https://GIS . stackexchange . com/questions/200679/convert-spatialprogonsdf-boundaries-to-spatial linesdf-keying-information-on-po

sp对象和栅格数据介绍(最完整)

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

合并多个SPDF

https://GIS . stackexchange . com/questions/155328/merging-multi-spatialpylondataforms-into-1-spdf-in-r

克里金插值

https://rpubs.com/nabilabd/118172

IDW插值

https://nceas . github . io/OSS-leaves/space-data-GIS-law/4-tues-space-analysis-in-r . html

Akima插值和ggplot2绘图

https://stack overflow . com/questions/50533738/best-of-method-of-spatial-interpolation-for-geographic-heat-等高线图

IDW插值和gplot2

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

Kknn插值和ggplot2

https://timogroshenbacher . ch/2018/03/classic-space-interpolation-with-r/

反向距离加权

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

将栅格图层对象转换为空间栅格/空间像素对象

https://GIS . stackexchange . com/questions/43707/how-to-product-spatial-grid-from-raster

KNN插值原理

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

美好的过去:

R_空间插值_必知必 会（一） ggplot2图集汇总（一） R_ggplot2地理信息可视化_史上最全（一） R语言中文社区2018年终文章整理（作者篇） R语言中文社区2018年终文章整理（类型篇）

微信官方账号可以在后台回复关键词学习

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