r语言插补后的数据 r语言空间插值

转载

数据大侠客 2024-05-23 13:55:30

文章标签 r语言插补后的数据 r生成新的dataframe 数据插值栅格 文章分类 R语言后端开发

作者：李誉辉

四川大学在读研究生

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 "E:/R_input_output/data_input/NanHai_9lines.csv", 
 9                       header = T, sep = ",")  
10colnames(Nine_lines) "long", "lat", "ID") 
11
12# 省会数据
13provinces "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 FALSE)
20
21x1 22xs1 23Chinaboundary_df 24
25Chinaboundary_df "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 FALSE)
31
32x2 33xs2 34Chinaprovinces_df 35Chinaprovinces_df "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 "+proj=aea +lat_1=25 +lat_2=50 +lon_0=105")
 8proj4string(Chinaboundary_sp)  9proj4string(Chinaprovinces_sp) 10
11# 将省会数据转化为sp对象
12provinces_sp 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 -3]
21  Sl_i 22  S_i list(Sl_i), ID = as.character(i))
23  list_lines[i] 24}
25
26Nine_lines_sp 27proj4string(Nine_lines_sp) 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 ",") # 
 7class(TEM_data1)
 8head(TEM_data1)
 9
10TEM_data1 1:10]  
11colnames(TEM_data1) "Station_Id", "latitude", "longitude", "altitude", # 重命名列
12                         "year", "month", "day", 
13                         "average_TEM", "highest_TEM", "lowest_TEM")  
14TEM_data2 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 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 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)23tm_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 7 Sr1 8 Srs1 list(Sr1), ID = "1") 9 SpP list(Srs1), 1:1)10 partmap_sp 11 Sr = SpP,12 data = data.frame(Names = "coords", row.names = row.names(SpP)))13 proj4string(partmap_sp) # 设定地图投影,sp转化为df后CRS丢失14 return(partmap_sp)15}1617# 分离台湾地图18Taiwan_sp 3.171, CRSchar = MyCRS) 19Taiwan_df 20Chinaboundary_noTaiwan_df 954.943, 2.903)) # 移除台湾21Sr1 22Srs1 list(Sr1), ID = "1") 23SpP list(Srs1), 1:1)24Chinaboundary_noTaiwan_sp 25 Sr = SpP,26 data = data.frame(Names = "coords", row.names = row.names(SpP)))27proj4string(Chinaboundary_noTaiwan_sp) 2829rm(Chinaboundary_noTaiwan_df, Sr1, Srs1, SpP) # 移除中途变量3031# 更改边框范围32TEM_sp@bbox

轮廓图：

1library(tmap)23tm_shape(shp = Chinaboundary_noTaiwan_sp) + 4 tm_borders(col = "magenta", lwd = 0.5) 56tm_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) "X", "Y")
10  coordinates(grd) "X", "Y")
11  gridded(grd) TRUE
12  fullgrid(grd) TRUE
13  proj4string(grd) 14  return(grd)
15}
16
17## 调用函数
18grd_TEM 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 1, # 反距离加权插值
 7                      locations = TEM_sp,
 8                      newdata = grd_TEM,
 9                      idp = 2) # 幂次为2
10TEM_raster # 栅格化
11TEM_mask # 筛选在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}1819tmap_TEM() 20

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。