问题描述
我想对德国各县的 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)
例如，我想评估当前的 PM10 值，直到 2020 年 6 月 1 日，我已经收到了来自德国联邦环境局的数据。数据如下： 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
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,Kreis,district,01001,DE.BB.BH
max values : DEU,DEU.9_1,Thüringen,DEU.9.9_1,Zwickau,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
创建 SpaceTime 对象：
month = 01:06
time = as.POSIXct(paste(month,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,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,])
> 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,0
variables : 14
names : GID_0,HASC_2,CC_2.2
min values : DEU,Alb-Donau-Kreis,Landkreis,08115,DE.BW.AD,8115
max values : DEU,DEU.1.6_1,Bodenseekreis,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：Plot STIDF
但是正如您所看到的，它无法正常工作。所以，我猜，我可能必须按月和县 ID 聚合：
pm10.f<-aggregate(PM10_f$TMW10,by = list(PM10_f$month,PM10_f$KRS),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,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,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 行和 polygon df 中的 6 个特征也遇到了同样的问题：
pm10.stf = STIDF(de,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,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 命令似乎只是从 polygon 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]),])
结果如下：
