R_空间插值_必知必会(一)

2019 年 2 月 18 日 R语言中文社区

作者:李誉辉  

四川大学在读研究生


1、简介

在气象等领域,空间插值非常重要,将观测站获取的数据汇总成点数据,
然后通过插值将点数据插值为栅格数据,再用地图boundary筛选出在boundary内的栅格。
最后将栅格数据添加到地图上。
本次教程会涉及到很多spsf的知识,十分详细
地图绘制采用2018年的新包tmap
,相比ggplot2更加简单快捷,
其语法类似,也是用
+号叠加图层,但是不能叠加其他图表。
同时也会增加
ggplot2中的绘制流程。
感谢南京信息工程学院大气物理学院的朋友提供的气象数据。
感谢张杰大佬在绘图过程中提供的帮助


2、中国地图数据

2.1

读取地图数据

分别读取:

  • 中国边界线地图数据

  • 中国省级地图数据

  • 南海九段线地图数据

  • 省会城市地图数据


 1rm(list = ls()); gc() # 清空内存
2library(sf)
3library(dplyr)
4library(magrittr)
5library(ggplot2)
6
7# 南海九段线数据
8Nine_lines <- read.csv(file = "E:/R_input_output/data_input/NanHai_9lines.csv"
9                       header = T, sep = ",")  
10colnames(Nine_lines) <- c("long""lat""ID"
11
12# 省会数据
13provinces <- read.csv(file = "E:/R_input_output/data_input/prov_centroids.csv"
14                           header = T)
15head(provinces)
16
17# 中国地图边界线官方数据
18path1 <- "E:/R_input_output/data_input/全国地图shp文件/1/bou1_4p.shp"
19Chinaboundary_sp <- rgdal::readOGR(dsn = path1, stringsAsFactors = FALSE)
20
21x1 <- Chinaboundary_sp@data
22xs1 <- data.frame(id = row.names(x1), x1)
23Chinaboundary_df <- fortify(Chinaboundary_sp)
24
25Chinaboundary_df <- full_join(Chinaboundary_df, xs1, by = "id")
26Chinaboundary_df %<>% filter(is.na(long) == FALSE & is.na(lat) == FALSE)# 去除空行
27
28# 中国省级地图官方数据
29path2 <- "E:/R_input_output/data_input/全国地图shp文件/1/bou2_4p.shp"
30Chinaprovinces_sp <- rgdal::readOGR(dsn = path2, stringsAsFactors = FALSE)
31
32x2 <- Chinaprovinces_sp@data
33xs2 <- data.frame(id = row.names(x2), x2)
34Chinaprovinces_df <- fortify(Chinaprovinces_sp)
35Chinaprovinces_df <- full_join(Chinaprovinces_df, xs2, by = "id")
36
37rm(path1, path2, x1, x2, xs1, xs2) # 移除中途变量


##           used (Mbgc trigger (Mbmax used (Mb)
## Ncells  519908 27.8    1180880 63.1   609142 32.6
## Vcells 1017814  7.8    8388608 64.0  1596878 12.2



## OGR data source with driver: ESRI Shapefile 
## Source: "E:\R_input_output\data_input\全国地图shp文件\1\bou1_4p.shp", layer: "bou1_4p"
## with 894 features
## It has 5 fields
## Integer64 fields read as strings:  BOU1_4M_ BOU1_4M_ID 
## OGR data source with driver: ESRI Shapefile 
## Source: "E:\R_input_output\data_input\全国地图shp文件\1\bou2_4p.shp", layer: "bou2_4p"
## with 925 features
## It has 7 fields
## Integer64 fields read as strings:  BOU2_4M_ BOU2_4M_ID


2.2

将dataframe转化为sp对象

 1library(sp)
2library(raster)
3library(magrittr)
4library(dplyr)
5
6# 统一坐标系
7MyCRS <- CRS("+proj=aea +lat_1=25 +lat_2=50 +lon_0=105")
8proj4string(Chinaboundary_sp) <- MyCRS
9proj4string(Chinaprovinces_sp) <- MyCRS
10
11# 将省会数据转化为sp对象
12provinces_sp <- SpatialPointsDataFrame(coords = cbind(x = provinces$x, y = provinces$y),
13                                 data = dplyr::select(provinces, name, shortname),
14                                 proj4string = MyCRS)
15
16# 将南海九段线数据转化为sp对象
17## SpatialLines()方法:
18list_lines <- list(NA)
19for (i in 1:length(unique(Nine_lines$ID))) {
20  l_i <- filter(Nine_lines, ID == i)[, -3]
21  Sl_i <- Line(l_i)
22  S_i <- Lines(list(Sl_i), ID = as.character(i))
23  list_lines[i] <- S_i
24}
25
26Nine_lines_sp <- SpatialLines(list_lines)
27proj4string(Nine_lines_sp) <- MyCRS
28rm(l_i, Sl_i, S_i, list_lines)


3、气象数据

3.1

读取气象数据

 1library(magrittr)
2library(dplyr)
3
4# 读取本地气象数据 
5path_TEM <- "E:/R_input_output/data_input/climate_data/2017_excel_processed/SURF_CLI_CHN_MUL_DAY-TEM-12001-201701.CSV"
6TEM_data1 <- read.csv(file = path_TEM, header = F,sep = ","
7class(TEM_data1)
8head(TEM_data1)
9
10TEM_data1 <- TEM_data1[,1:10]  
11colnames(TEM_data1) <- c("Station_Id""latitude""longitude""altitude"# 重命名列
12                         "year""month""day"
13                         "average_TEM""highest_TEM""lowest_TEM")  
14TEM_data2 <- transmute(TEM_data1, 
15                        lat = latitude / 100
16                        long = longitude / 100
17                        alt = altitude / 10
18                        aver_TEM = average_TEM / 10
19                        high_TEM = highest_TEM / 10
20                        low_TEM = lowest_TEM / 10)
21
22TEM_data3 <-  cbind(dplyr::select(TEM_data1, "Station_Id""year""month""day"),
23                               TEM_data2)
24
25head(TEM_data3)
26nrow(TEM_data3)
27
28length(unique(TEM_data3$month)) # 只有1个元素,都是1月的数据
29unique(TEM_data3$year) # 只有1个元素,都是2017年的
30length(unique(TEM_data3$day)) # 31个元素,刚好31天
31TEM_1th <- subset(TEM_data3, day == 1# 筛选1月1日的数据
32
33rm(path_TEM, TEM_data1, TEM_data2, TEM_data3) # 移除数据生成过程中的变量


## [1] "data.frame"



## [1] 26009
## [1] 1
## [1] 2017
## [1] 31


3.2

将dataframe转化为sp对象

 1library(sp)
2library(raster)
3library(magrittr)
4library(dplyr)
5
6# 统一坐标系
7MyCRS <- CRS("+proj=aea +lat_1=25 +lat_2=50 +lon_0=105")
8# 将温度数据转化为sp对象
9TEM_sp <- SpatialPointsDataFrame(coords = cbind(x = TEM_1th$long, y = TEM_1th$lat),
10                                 data = dplyr::select(TEM_1th, aver_TEM),
11                                 proj4string = MyCRS)
12


4、空间插值边走边学

4.1

散点地图

根据散点图确定栅格边界范围。
从散点图中,可以看出,没有台湾地区的数据,所以将台湾地图分离出来。

1library(tmap)
2
3tm_shape(shp = Chinaboundary_sp) + 
4  tm_polygons(col = "yellow", border.col = "cyan") + 
5  tm_borders(lwd = 0.5) + 
6  tm_shape(shp = TEM_sp) + 
7  tm_dots(col = "magenta", size = 0.2, shape = 21



4.2

分离台湾地图

很多情况下,中国官方发布的数据中并不包含台湾,
但地图需要有台湾的轮廓线,
所以需要将台湾部分分离出来
分离原理:
大陆主体部分
AREA == 954.943,
海南岛:
AREA == 2.903,
台湾:
AREA == 954.943

1library(sp)
2library(dplyr)
3
4# 生成一个函数,将dataframe转化为sp对象,
5map_separation <- function(map_df, area, CRSchar) # x,y指定经纬度
6  map_subset <-  subset(map_df, AREA == area)
7  Sr1 <- Polygon(cbind(map_subset$long, map_subset$lat)) 
8  Srs1 <- Polygons(list(Sr1), ID = "1"
9  SpP <- SpatialPolygons(Srl = list(Srs1), 1:1)
10  partmap_sp <- SpatialPolygonsDataFrame(
11    Sr = SpP,
12    data = data.frame(Names = "coords", row.names = row.names(SpP)))
13  proj4string(partmap_sp) <- CRSchar # 设定地图投影,sp转化为df后CRS丢失
14  return(partmap_sp)
15}
16
17# 分离台湾地图
18Taiwan_sp <- map_separation(Chinaboundary_df, 3.171, CRSchar = MyCRS) 
19Taiwan_df <- fortify(Taiwan_sp) 
20Chinaboundary_noTaiwan_df <-  subset(Chinaboundary_df, AREA %in% c(954.9432.903)) # 移除台湾
21Sr1 <- Polygon(cbind(Chinaboundary_noTaiwan_df$long, Chinaboundary_noTaiwan_df$lat)) 
22Srs1 <- Polygons(list(Sr1), ID = "1"
23SpP <- SpatialPolygons(Srl = list(Srs1), 1:1)
24Chinaboundary_noTaiwan_sp <- SpatialPolygonsDataFrame(
25  Sr = SpP,
26  data = data.frame(Names = "coords", row.names = row.names(SpP)))
27proj4string(Chinaboundary_noTaiwan_sp) <- MyCRS
28
29rm(Chinaboundary_noTaiwan_df, Sr1, Srs1, SpP) # 移除中途变量
30
31# 更改边框范围
32TEM_sp@bbox <- Chinaboundary_sp@bbox

轮廓图:

1library(tmap)
2
3tm_shape(shp = Chinaboundary_noTaiwan_sp) + 
4  tm_borders(col = "magenta", lwd = 0.5
5
6tm_shape(shp = Taiwan_sp) + 
7  tm_borders(col = "magenta", lwd = 0.5



4.3

生成空栅格

 1ibrary(gstat) # 内含插值函数
2library(sp)    # 内含spsample函数
3library(raster)
4
5# 创建空栅格,n表示栅格数量
6## 首先生成空栅格函数
7grd_empty <- function(wheather_sp, n) {
8  grd <- as.data.frame(spsample(wheather_sp, "regular", n=n))
9  names(grd) <- c("X""Y")
10  coordinates(grd) <- c("X""Y")
11  gridded(grd) <- TRUE
12  fullgrid(grd) <- TRUE
13  proj4string(grd) <- proj4string(wheather_sp)
14  return(grd)
15}
16
17## 调用函数
18grd_TEM <- grd_empty(wheather_sp = TEM_sp, n = 4000000# 栅格总数量


4.4

IDW(反距离加权)插值

通过插值给栅格赋值。
反距离加权插值,是近期做大数据显示时使用的插值方法,很好用的插值方法。
距离cell约近的数据点,对其影响越大。需要指定幂值P,

P默认为2,一般为0.5到3均可获得合理的结果。
通过定义更高的幂值,可进一步强调最近点。
因此,邻近数据将受到更大影响,表面会变得更加详细(更不平滑)。
随着幂数的增大,内插值将逐渐接近最近采样点的值。
指定较小的幂值将对距离较远的周围点产生更大的影响,从而导致平面更加平滑。

 1library(gstat) # 内含插值函数
2library(sp)    # 内含spsample函数
3library(raster)
4
5# idw插值(指定幂次为2, 即idp = 2)创建raster
6TEM_idw <- gstat::idw(formula = aver_TEM ~ 1# 反距离加权插值
7                      locations = TEM_sp,
8                      newdata = grd_TEM,
9                      idp = 2# 幂次为2
10TEM_raster <- raster(TEM_idw) # 栅格化
11TEM_mask <- mask(TEM_raster, Chinaboundary_noTaiwan_sp) # 筛选在boundary范围内的栅格
12rm(TEM_idw, TEM_raster)


## [inverse distance weighted interpolation]


5、tmap画图
  • tm_shape()用于传入数据,支持spsf对象。不支持dataframe对象。

  • tm_raster()用于添加栅格图层。

  • tm_border()用于添加边界线。

  • tm_lines()用于添加多个线段。

  • tm_polygons()用于添加多边形。

  • tm_text()用于添加文本。

  • tm_legend()用于图例格式设定。

tmap支持RColorBrewer中所有色板。色板名称前加-就能颜色标度反向。
总之与
ggplot2类似,若对ggplot2画地图不熟悉的,可以参考R_ggplot2地理信息可视化_史上最全
ggplot2其它图形不属性的可以参考R_ggplot2基础(四),该教程是一个连载,在下一篇文末有其它章节链接。

 1library(tmap)
2
3tmap_TEM <- function(temperature_data = TEM_mask) {
4  tm_shape(shp = temperature_data) + 
5    tm_raster(n=10, palette = "-PuOr"legend.reverse = TRUE, # 负色板
6              auto.palette.mapping = FALSE,
7              title="中国平均气温\n2017年1月1日 \n(单位:摄氏度)") + 
8    tm_shape(shp = Chinaprovinces_sp) +  # 省级地图数据
9    tm_borders(col = "black"lwd = 0.5) + 
10    tm_shape(shp = Nine_lines_sp) +   # 南海九段线数据
11    tm_lines(col = "blue"lwd = 3) + 
12    tm_shape(shp = Taiwan_sp) +  # 给台湾地图改变颜色,避免与业务数据颜色相同
13    tm_polygons(col = "grey") + 
14    tm_shape(shp = provinces_sp) +  # 省会名称数据
15    tm_text(text = "shortname"col = "green"size = 0.5) + # 注意变量名加引号
16    tm_legend(legend.outside = TRUE)
17}
18
19tmap_TEM() 
20


往期精彩:


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

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

经过我的不懈努力,体重终于掉了两斤↓

其中1斤是头发↓

登录查看更多
0

相关内容

【ICML2020】对比多视角表示学习
专知会员服务
52+阅读 · 2020年6月28日
【MIT-ICML2020】图神经网络的泛化与表示的局限
专知会员服务
41+阅读 · 2020年6月23日
专知会员服务
59+阅读 · 2020年3月19日
【ICLR-2020】网络反卷积,NETWORK DECONVOLUTION
专知会员服务
37+阅读 · 2020年2月21日
CVPR 2019 | 34篇 CVPR 2019 论文实现代码
AI科技评论
21+阅读 · 2019年6月23日
R语言机器学习:xgboost的使用及其模型解释
R语言中文社区
10+阅读 · 2019年5月6日
卷积神经网络四种卷积类型
炼数成金订阅号
18+阅读 · 2019年4月16日
超详细干货 | 三维语义分割概述及总结
计算机视觉life
32+阅读 · 2019年3月19日
一文带你读懂 DeconvNet 上采样层(语义分割)
AI研习社
26+阅读 · 2019年3月16日
基于R语言进行Box-Cox变换
R语言中文社区
45+阅读 · 2018年11月19日
基于混合张量分解提升扩张卷积网络
论智
11+阅读 · 2018年2月11日
Capsule Networks解析
机器学习研究会
10+阅读 · 2017年11月12日
Learning Dynamic Routing for Semantic Segmentation
Arxiv
8+阅读 · 2020年3月23日
Learning Implicit Fields for Generative Shape Modeling
Arxiv
10+阅读 · 2018年12月6日
Arxiv
8+阅读 · 2018年5月15日
VIP会员
相关资讯
CVPR 2019 | 34篇 CVPR 2019 论文实现代码
AI科技评论
21+阅读 · 2019年6月23日
R语言机器学习:xgboost的使用及其模型解释
R语言中文社区
10+阅读 · 2019年5月6日
卷积神经网络四种卷积类型
炼数成金订阅号
18+阅读 · 2019年4月16日
超详细干货 | 三维语义分割概述及总结
计算机视觉life
32+阅读 · 2019年3月19日
一文带你读懂 DeconvNet 上采样层(语义分割)
AI研习社
26+阅读 · 2019年3月16日
基于R语言进行Box-Cox变换
R语言中文社区
45+阅读 · 2018年11月19日
基于混合张量分解提升扩张卷积网络
论智
11+阅读 · 2018年2月11日
Capsule Networks解析
机器学习研究会
10+阅读 · 2017年11月12日
Top
微信扫码咨询专知VIP会员