r - 时空包的问题
问题描述
我想对德国各县的 PM10 进行每月时空分析并绘制它们。稍后我想分析不同的回归模型。但是我无法创建一个时空物体，我需要进一步分析和其他我将要研究的问题。所以，我首先开始了解方法和包，尽我所能，我被困在我无法创建适当的时空对象的点上。
我以以下可重现代码为指导（来源：https ://edzer.github.io/UseR2016/ ）：
data("Produc", package = "plm")
Produc[1:5,1:9]
library(maps)
states.m = map('state', plot=FALSE, fill=TRUE)
IDs <- sapply(strsplit(states.m$names, ":"), function(x) x[1])
library(maptools)
states = map2SpatialPolygons(states.m, IDs=IDs)
yrs = 1970:1986
time = as.POSIXct(paste(yrs, "-01-01", sep=""), tz = "GMT")
time
library(spacetime)
Produc.st = STFDF(states[-8], time, Produc[order(Produc[2], Produc[1]),])
library(RColorBrewer)
stplot(Produc.st[,,"unemp"], yrs, col.regions = brewer.pal(9, "YlOrRd"), cuts = 9)
例如，我想评估截至 2020 年 6 月 1 日的县级当前 PM10 值，为此我收到了德国联邦环境署的数据。数据如下： PM10 是我的df，感兴趣的值是TMW，这是 PM10 的每日平均测量值。
PM10[sample(nrow(PM10),10),]
# A tibble: 10 x 9
Station Komponente Datum TYPEOFAREA TYPEOFSTATION TMW TMW_R TypeOfData Lieferung
<chr> <chr> <date> <chr> <chr> <dbl> <dbl> <chr> <chr>
1 DENI051 PM10 2020-02-28 ländliches Gebiet Hintergrund 5.40 5 S M
2 DETH095 PM10 2020-05-12 städtisches Gebiet Hintergrund 9.74 10 S M
3 DEBY118 PM10 2020-04-30 städtisches Gebiet Hintergrund 5.27 5 S M
4 DEBY072 PM10 2020-05-03 ländlich regional Hintergrund 8.43 8 S M
5 DEHE060 PM10 2020-06-01 ländlich regional Hintergrund 9.43 9 S M
6 DEBW087 PM10 2020-05-28 ländlich regional Hintergrund 11.0 11 S M
7 DEBW038 PM10 2020-03-11 städtisches Gebiet Hintergrund 4.28 4 S M
8 DENW065 PM10 2020-01-10 ländlich regional Hintergrund 2.16 2 S M
9 DENW096 PM10 2020-05-17 vorstädtisches Gebiet Hintergrund 13.2 13 T M
10 DEHE050 PM10 2020-04-20 ländliches Gebiet Hintergrund 8.20 8 S D
然后我从https://gadm.org/download_country_v3.html --> Germany --> R(sp) --> level2下载了一个 sp 文件
其中包含德国县级地图，如下所示：
> de
class : SpatialPolygonsDataFrame
features : 403
extent : 5.866251, 15.04181, 47.27012, 55.05653 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 13
names : GID_0, NAME_0, GID_1, NAME_1, NL_NAME_1, GID_2, NAME_2, VARNAME_2, NL_NAME_2, TYPE_2, ENGTYPE_2, CC_2, HASC_2
min values : DEU, Germany, DEU.1_1, Baden-Württemberg, NA, DEU.1.1_1, Ahrweiler, NA, NA, Kreis, District, 01001, DE.BB.BH
max values : DEU, Germany, DEU.9_1, Thüringen, NA, DEU.9.9_1, Zwickau, NA, NA, Water body, Water body, 16077, DE.TH.WR
由于我的 df 不包括县级的地理配准而是站代码，因此我已将此信息添加到数据集中。我的 sp 文件中的县 ID 是CC_2，如果 ID 有四位数，它是一个以 0 开头的五位数代码。例子：
de$CC_2
[1] "08425" "08211" "08426" "08115" "12065" "12066" "12067"
我猜的第一个问题是，当我通过站代码将地理信息添加到我的df时，我在 df 中得到了我的CC_2，如下所示：
> PM10_m[sample(nrow(PM10_m),3),]
Station Komponente Datum TYPEOFAREA TYPEOFSTATION TMW TMW_R TypeOfData Lieferung CC_2
11448 DEBW081 PM10 2020-06-07 städtisches Gebiet Hintergrund 6.775362 7 T M 8212
1566 DEBB066 PM10 2020-04-19 ländlich regional Hintergrund 11.162500 11 S M 12061
7174 DEBW027 PM10 2020-03-20 städtisches Gebiet Hintergrund 34.791667 35 S M 8415
如您所见，缺少了四位 ID 开头的 0，因此我检查了变量的结构：
str(PM10_m$CC_2)
chr [1:47350] "12062" "12062" "12062" "12062" "12062" "12062" "12062" "12062" "12062" "12062" "12062" "12062" "12062" ...
str(de$CC_2)
chr [1:403] "08425" "08211" "08426" "08115" NA "08435" "08315" "08235" "08316" "08236" "08116" "08311" "08237" "08117" ...
因此，两者都是chr，但如果将它们匹配每四位数字 ID 将不匹配！所以，我过去常常通过将两个变量都设为数字来处理这个问题。在这一点上，我不确定我这样做是否正确。
> PM10_m$CC_2<-as.numeric(PM10_m$CC_2)
> de$CC_2.2<-as.numeric(de$CC_2)
在合并它们之前，我曾经按县 ID 和日期聚合 PM10_m df。
PM10_aggr<-aggregate(PM10_m$TMW, by = list(PM10_m$Datum, PM10_m$CC_2), FUN="mean", na.rm=T)
我现在合并了 df 和多边形 df de，看看它是否有效。
de_t<-merge(de,PM10_aggr, by.x="CC_2.2", by.y="CC_2", na.rm=T,duplicateGeoms=TRUE)
据我所见，它匹配正确： Plotting with tmap
现在，我开始创建一个时空对象，按照指南中的步骤（见开头）：
首先，我将月份添加到我的 df PM10_aggr
PM10_f<-PM10_aggr
PM10_f$month<-strftime(PM10_aggr$date, format = "%m")
> PM10_f[sample(nrow(PM10_f),4),]
date CC_2 TMW10 month
26303 2020-04-04 13062 6.136208 04
24703 2020-05-12 12072 7.506250 05
4808 2020-03-16 3452 13.933222 03
30502 2020-04-17 16051 30.121002 04
创建时空对象：
month = 01:06
time = as.POSIXct(paste(month, "-01-01", sep=""), tz = "GMT")
time
[1] "0001-01-01 GMT" "0002-01-01 GMT" "0003-01-01 GMT" "0004-01-01 GMT" "0005-01-01 GMT" "0006-01-01 GMT"
它不像指南中那样工作，但据我所知，它只是创建和分类时间对象。所以，我走在指南前面：
library(spacetime)
pm10.st = STFDF(de, time, PM10_f[order(PM10_f[4], PM10_f[1]),])
Error in validityMethod(object) :
nrow(object@data) == length(object@sp) * nrow(object@time) is not TRUE
我读到命令STFDF不能处理丢失的地理点，我必须改用命令STIDF。
所以，这就是我得到的：
pm10.st = STIDF(de, time, PM10_f[order(PM10_f[4], PM10_f[1]),])
> pm10.st
An object of class "STIDF"
Slot "data":
date KRS TMW10 month month1
1 2020-01-01 1002 33.34608 01 1
183 2020-01-01 1003 81.06596 01 1
365 2020-01-01 1051 53.14400 01 1
547 2020-01-01 1053 34.36517 01 1
729 2020-01-01 1054 NaN 01 1
911 2020-01-01 1057 32.04604 01 1
Slot "sp":
class : SpatialPolygonsDataFrame
features : 6
extent : 8.108812, 10.24141, 47.5024, 48.86768 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 14
names : GID_0, NAME_0, GID_1, NAME_1, NL_NAME_1, GID_2, NAME_2, VARNAME_2, NL_NAME_2, TYPE_2, ENGTYPE_2, CC_2, HASC_2, CC_2.2
min values : DEU, Germany, DEU.1_1, Baden-Württemberg, NA, DEU.1.1_1, Alb-Donau-Kreis, NA, NA, Landkreis, District, 08115, DE.BW.AD, 8115
max values : DEU, Germany, DEU.1_1, Baden-Württemberg, NA, DEU.1.6_1, Bodenseekreis, NA, NA, Water body, Water body, 08435, DE.BW.BR, 8435
Slot "time":
timeIndex
0001-01-01 1
0002-01-01 2
0003-01-01 3
0004-01-01 4
0005-01-01 5
0006-01-01 6
Slot "endTime":
[1] "0001-01-01 GMT" "0002-01-01 GMT" "0003-01-01 GMT" "0004-01-01 GMT" "0005-01-01 GMT" "0006-01-01 GMT"
当我看到该命令仅从 df 中取出 6 行并且仅与多边形 df 的 6 个特征匹配时，我真的很惊讶。我可以绘制这个 STIDF：绘制 STIDF
但正如你所看到的，它不能正常工作。所以，我猜，我可能必须按月份和县 ID 汇总：
pm10.f<-aggregate(PM10_f$TMW10, by = list(PM10_f$month, PM10_f$KRS),FUN="mean", na.rm=T)
> str(pm10.f)
'data.frame': 1092 obs. of 3 variables:
$ month: chr "01" "02" "03" "04" ...
$ CID : num 1002 1002 1002 1002 1002 ...
$ MMW10: num 13.3 11.1 14.2 16.1 12.4 ...
### CID is the County ID ###
> pm10.f[sample(nrow(pm10.f),5),]
month CID MMW10
234 06 5158 16.637490
704 02 9775 11.083747
1030 04 16055 18.934881
842 02 13054 8.594628
513 03 8121 16.9119
所以，我再次尝试使用STIDF 命令：
pm10.stf = STIDF(de, time, pm10.f[order(pm10.f[1], pm10.f[1]),])
> pm10.stf
An object of class "STIDF"
Slot "data":
month CID MMW10
1 01 1002 13.31264
7 01 1003 17.81540
13 01 1051 17.67919
19 01 1053 12.99228
25 01 1054 NaN
31 01 1057 14.71878
Slot "sp":
class : SpatialPolygonsDataFrame
features : 6
extent : 8.108812, 10.24141, 47.5024, 48.86768 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 14
names : GID_0, NAME_0, GID_1, NAME_1, NL_NAME_1, GID_2, NAME_2, VARNAME_2, NL_NAME_2, TYPE_2, ENGTYPE_2, CC_2, HASC_2, CC_2.2
min values : DEU, Germany, DEU.1_1, Baden-Württemberg, NA, DEU.1.1_1, Alb-Donau-Kreis, NA, NA, Landkreis, District, 08115, DE.BW.AD, 8115
max values : DEU, Germany, DEU.1_1, Baden-Württemberg, NA, DEU.1.6_1, Bodenseekreis, NA, NA, Water body, Water body, 08435, DE.BW.BR, 8435
Slot "time":
timeIndex
0001-01-01 1
0002-01-01 2
0003-01-01 3
0004-01-01 4
0005-01-01 5
0006-01-01 6
Slot "endTime":
[1] "0001-01-01 GMT" "0002-01-01 GMT" "0003-01-01 GMT" "0004-01-01 GMT" "0005-01-01 GMT" "0006-01-01 GMT"
我遇到了同样的问题，再次只有 6 个随机行与 6 个县匹配：plot STIDF 2
即使我删除了order 命令，我也遇到了同样的问题，只有 df 的 6 行和多边形 df的6 个特征：
pm10.stf = STIDF(de, time, pm10.f)
> pm10.stf
An object of class "STIDF"
Slot "data":
month CID MMW10
1 01 1002 13.31264
2 02 1002 11.10590
3 03 1002 14.19649
4 04 1002 16.10512
5 05 1002 12.38511
6 06 1002 13.10104
Slot "sp":
class : SpatialPolygonsDataFrame
features : 6
extent : 8.108812, 10.24141, 47.5024, 48.86768 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 14
names : GID_0, NAME_0, GID_1, NAME_1, NL_NAME_1, GID_2, NAME_2, VARNAME_2, NL_NAME_2, TYPE_2, ENGTYPE_2, CC_2, HASC_2, CC_2.2
min values : DEU, Germany, DEU.1_1, Baden-Württemberg, NA, DEU.1.1_1, Alb-Donau-Kreis, NA, NA, Landkreis, District, 08115, DE.BW.AD, 8115
max values : DEU, Germany, DEU.1_1, Baden-Württemberg, NA, DEU.1.6_1, Bodenseekreis, NA, NA, Water body, Water body, 08435, DE.BW.BR, 8435
Slot "time":
timeIndex
0001-01-01 1
0002-01-01 2
0003-01-01 3
0004-01-01 4
0005-01-01 5
0006-01-01 6
Slot "endTime":
[1] "0001-01-01 GMT" "0002-01-01 GMT" "0003-01-01 GMT" "0004-01-01 GMT" "0005-01-01 GMT" "0006-01-01 GMT"
我在df中有一个县的 6 行，但有不同的 6 个多边形特征。似乎STIDF命令只是从多边形 df中获取前 6 个多边形。
解决方案
首先，我注意到我的 shapefile 包含的元素比实际的区域多。这是因为 shapefile 包含“DoubleGeoms”。所以我将shapefile聚合如下：
raster::aggregate(de, by="AGS")
然后我突然想到我的想法有逻辑错误。所以我有 401 个区，实际上有 6 个测量时间（6 个月），所以我的数据框应该有 401*6=2406 行。这意味着我必须调整我的数据框。所以我把401个区都扩大了：
df<-tidyr::expand_grid(KRS=df$KRS,1:6)
使用“merge”命令按地区和月份将变量添加到新数据帧后，我现在可以使用“ spacetime ”包中的“ STFDF ”命令：
df.stf <- STFDF(de2, time, df[order(df[2], df[1]),])
这是结果：
推荐阅读
- firebase - 新手bigquery如何从firebase事件中选择多个参数（模式记录，重复）
- perl - 使用 Date::Calc 的年份和日期
- javascript - 我怎样才能让统计框显示在我的threejs画布上？
- swift - 将 viewDidLoad() 或 ViewDidAppear() 中的变量传递给按钮操作
- php - 谷歌地图 api 标记、sql 和 xml
- rest - Servicenow 与 SAP 的集成
- c - 为什么我在 C 中遇到分段错误？
- excel - Excel VBA宏读取具有不同文本的一列
- odata - sap.m.table 未更新
- google-bigquery - BigQuery 是否可以使用 aes_decrypt 函数？