#r #spatial #sf #sp #rgdal
#r #пространственный #sf #sp #rgdal
Вопрос:
Я хочу сделать ежемесячный пространственно-временной анализ 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-06-01 ежемесячно на уровне округа для этого я получил данные от Федерального агентства по охране окружающей среды Германии. Данные выглядят следующим образом:
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
Затем я загрузил sp-файл из https://gadm.org/download_country_v3.html -> Германия -> R(sp) -> уровень 2
который содержит карту Германии на уровне округа и выглядит следующим образом:
> 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-файле — CC_2, который представляет собой пятизначный код, начинающийся с 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
As you can see, the 0 in the beginning of a four digit ID is missing, so I checked the structure of the variables:
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" ...
So, both are chr but if would match them up every four digit ID would not match! So, I used to handle this by making both variables as a numeric. At this point I am not sure, if it is right when I am doing it this way.
> PM10_m$CC_2<-as.numeric(PM10_m$CC_2)
> de$CC_2.2<-as.numeric(de$CC_2)
Before I merge them, I used to aggregate the PM10_m df by county ID and date.
PM10_aggr<-aggregate(PM10_m$TMW, by = list(PM10_m$Datum, PM10_m$CC_2), FUN="mean", na.rm=T)
I merged now the df and the polygon df de, to see if it worked.
de_t<-merge(de,PM10_aggr, by.x="CC_2.2", by.y="CC_2", na.rm=T,duplicateGeoms=TRUE)
As far I can see, it matched properly:
Plotting with tmap
Now, I started to create a spacetime object, following the steps like in the guide (see in the beginning):
First I added the month into my 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
Creating SpaceTime object:
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"
It worked not like in the guide but as far as I understand, it is just creating and categorial time object. So, I steped forward of the guide:
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
I read that the command STFDF can’t work with missing geopoints and that I have to use the command STIDF instead.
So, this is what I get:
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"
I was really suprised when I saw, that the command just took 6 rows from the df and matches with just 6 features of the polygon df. I can Plot this STIDF:Plot STIDF
But as you can see it does not worked properly. So, I guessed, may I have to aggregate by month and county 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
So, I tried again with the STIDF command:
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, у меня возникнут те же проблемы только с 6 строками из df и 6 объектами из polygon 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"
Я получил 6 строк одного округа в df, но разные 6 полигональных объектов. Кажется, что команда STIDF просто берет первые 6 полигонов из polygon df.
Комментарии:
1. Хорошее описание того, что вы сделали, но его трудно понять суть вопроса. Не могли бы вы изложить вопрос / желаемый результат в верхней части сообщения в виде краткого резюме, прежде чем описывать, что вы сделали?
2. @AndreWildberg я отредактировал свой вопрос и добавил свою цель, и я надеюсь, что это касается моих действий. Спасибо за ваш совет!
Ответ №1:
Во-первых, я заметил, что в моем шейп-файле больше элементов, чем фактических районов. Это связано с тем, что шейп-файл содержит «двойные геометрические фигуры». Итак, я объединил шейп-файл следующим образом:
raster::aggregate(de, by="AGS")
Затем мне пришло в голову, что у меня логическая ошибка в мышлении. Итак, у меня 401 округ и практически 6 раз измерения (6 месяцев), поэтому в моем фрейме данных должно быть 401 * 6 = 2406 строк. Это означает, что мне пришлось настроить свой фрейм данных. Итак, я взял 401 округ и расширил их:
df<-tidyr::expand_grid(KRS=df$KRS,1:6)
После добавления переменных в новый фрейм данных с помощью команды «объединить» по районам и месяцам, теперь я мог бы использовать команду «STFDF» из пакета «spacetime«:
df.stf <- STFDF(de2, time, df[order(df[2], df[1]),])
И это результат:
