欢迎关注天善智能，我们是专注于商业智能BI，人工智能AI，大数据分析与挖掘领域的垂直社区，学习，问答、求职一站式搞定！
对商业智能BI、大数据分析挖掘、机器学习，python，R等数据领域感兴趣的同学加微信：tstoutiao，邀请你进入数据爱好者交流群，数据爱好者们都在这儿。
作者：李誉辉
四川大学在读研究生
1、简介
在气象等领域，空间插值非常重要，将观测站获取的数据汇总成点数据，
然后通过插值将点数据插值为栅格数据，再用地图boundary筛选出在boundary内的栅格。
最后将栅格数据添加到地图上。
本次教程会涉及到很多sp和sf的知识，十分详细。
地图绘制采用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 (Mb) gc trigger (Mb) max 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.943, 2.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()用于传入数据，支持sp和sf对象。不支持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
往期精彩：
-
ggplot2图集汇总（一）
-
R_ggplot2地理信息可视化_史上最全（一）
-
R语言中文社区2018年终文章整理（作者篇）
-
R语言中文社区2018年终文章整理（类型篇）
公众号后台回复关键字即可学习
回复 爬虫 爬虫三大案例实战
回复 Python 1小时破冰入门
回复 数据挖掘 R语言入门及数据挖掘
回复 人工智能 三个月入门人工智能
回复 数据分析师 数据分析师成长之路
回复 机器学习 机器学习的商业应用
回复 数据科学 数据科学实战
回复 常用算法 常用数据挖掘算法
网友评论