时空包的问题
problems with spacetime package
我想对德国各县的 PM10 进行月度 space 时间分析并绘制它们。稍后我想分析不同的回归模型。但是我无法创建一个 spacetime 对象，我需要它来进一步分析和我将要处理的其他研究问题。所以，我首先开始了解方法和包，尽我所能，我被困在这一点上，我无法创建一个合适的 spacetime 对象。
我正在将以下可重现代码作为指南（来源：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 下载了一个 sp 文件 --> 德国 --> R(sp) --> level2
其中包含德国县级地图，如下所示：
> 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 时，我得到了 CC_2 在 df 中是这样的：
> 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
现在，我开始创建一个 spacetime 对象，按照 指南 中的步骤（见开头）：
首先我将月份添加到我的 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, "-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：Plot 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 行和来自 [=80] 的 6 个特征=]多边形 df:
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 命令似乎只是从 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)
在按地区和月份使用“合并”命令将变量添加到新数据帧后，我现在可以使用“”中的“STFDF”命令时空”包：
df.stf <- STFDF(de2, time, df[order(df[2], df[1]),])
这是结果：
我想对德国各县的 PM10 进行月度 space 时间分析并绘制它们。稍后我想分析不同的回归模型。但是我无法创建一个 spacetime 对象，我需要它来进一步分析和我将要处理的其他研究问题。所以，我首先开始了解方法和包，尽我所能，我被困在这一点上，我无法创建一个合适的 spacetime 对象。
我正在将以下可重现代码作为指南（来源：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 下载了一个 sp 文件 --> 德国 --> R(sp) --> level2
其中包含德国县级地图，如下所示：
> 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 时，我得到了 CC_2 在 df 中是这样的：
> 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
现在，我开始创建一个 spacetime 对象，按照 指南 中的步骤（见开头）：
首先我将月份添加到我的 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, "-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：Plot 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 行和来自 [=80] 的 6 个特征=]多边形 df:
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 命令似乎只是从 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)
在按地区和月份使用“合并”命令将变量添加到新数据帧后，我现在可以使用“”中的“STFDF”命令时空”包：
df.stf <- STFDF(de2, time, df[order(df[2], df[1]),])
这是结果：