|
作者：李誉辉 四川大学在读研究生
在气象等领域，空间插值非常重要，将观测站获取的数据汇总成点数据，
然后通过插值将点数据插值为栅格数据，再用地图boundary筛选出在boundary内的栅格。
最后将栅格数据添加到地图上。
本次教程会涉及到很多sp和sf的知识，十分详细。
地图绘制采用2018年的新包tmap，相比ggplot2更加简单快捷，
其语法类似，也是用+号叠加图层，但是不能叠加其他图表。
同时也会增加ggplot2中的绘制流程。
感谢南京信息工程学院大气物理学院的朋友提供的气象数据。
感谢张杰大佬在绘图过程中提供的帮助。
分别读取： 中国边界线地图数据
中国省级地图数据
南海九段线地图数据
省会城市地图数据
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
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)
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] 26009
## [1] 1
## [1] 2017
## [1] 31
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
根据散点图确定栅格边界范围。
从散点图中，可以看出，没有台湾地区的数据，所以将台湾地图分离出来。 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)
很多情况下，中国官方发布的数据中并不包含台湾，
但地图需要有台湾的轮廓线,
所以需要将台湾部分分离出来
分离原理：
大陆主体部分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)
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) # 栅格总数量
通过插值给栅格赋值。
反距离加权插值，是近期做大数据显示时使用的插值方法，很好用的插值方法。
距离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]
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
|
