In volkerschmid/luftdaten: Analysing data from luftdaten.info

Only Munich, Only a couple of hours

library(osmdata)
query <- opq(bbox='Munich, Germany')
bbox<-as.numeric(strsplit(query$bbox,",")[[1]])

library(RSQLite)
db<-dbConnect(SQLite(),"../db/CopyOf2018-SDS011.sqlite")
P1<-dbGetQuery(db,paste0("SELECT dataP.locid,time(dataP.timestamp),locid.lat,locid.lon,dataP.P1 from dataP,locid where date(timestamp)='2018-02-28' and locid.id=dataP.locid and locid.lat<",bbox[3]," AND locid.lat>",bbox[1]," AND locid.lon>",bbox[2], " AND locid.lon<",bbox[4]))
print(dim(P1))
locid<-dbGetQuery(db,paste0("SELECT locid.id,locid.lat,locid.lon from dataP,locid where date(timestamp)='2018-02-28' and locid.id=dataP.locid and locid.lat<",bbox[3]," AND locid.lat>",bbox[1]," AND locid.lon>",bbox[2], " AND locid.lon<",bbox[4]))
locid<-unique(locid)
dbDisconnect(db)

We have 114 locations with 62931 observations. Plotting locations:

if (0){
#Notice: Osmar with local file is very slow, because file os so big.
      ctown <- get_osm(bb, source = src)
      for (j in c("motorway","trunk","primary","secondary","tertiary"))
        {
        road <- find(ctown, way(tags(k == "highway" & v==j)))
        road <- find_down(ctown, way(road))
        road <- subset(ctown, ids = road)
        plot_ways(road, add=TRUE)
        }
        road <- find(ctown, way(tags(k == "highway" & v=="residential")))
        road <- find_down(ctown, way(road))
        road <- subset(ctown, ids = road)
        plot_ways(road, add=TRUE, col=grey(.6))
}

luftdaten::plot_loc(locid,range=as.numeric(bbox[c(2,4,1,3)]))
query <- add_osm_feature(query, key='highway', value='!residential')
test<-osmdata_sp(query)
sp::plot(test$osm_lines, add=TRUE)

(Reduce time, try afternoon 12pm top 4pm) compute time as POSIX

times<-as.POSIXct(P1$`time(dataP.timestamp)`,format="%H:%M:%S")
#P1<-P1[times<paste(Sys.Date(),"15:00:00")&times>paste(Sys.Date(),"12:00:00"),]
#times<-times[times<paste(Sys.Date(),"15:00:00")&times>paste(Sys.Date(),"12:00:00")]
print(dim(P1))
print(length(unique((P1$locid))))

construct data for spTimer, ten-minutes means (gives 6x24=144 time points)

locations<-sort(unique(locid$id))[-55]
L<-length(locations)
TI<-144
data<-rep(NA,L*TI)
counter<-0
# bad coding ahead
for (j in 1:L)
{
    find2<-P1$locid==locations[j]
      for (t in 1:T)
      {
        find<-(times<as.POSIXct(paste(Sys.Date(),"00:00:00"))+600*t)
        find<-find&(times>=as.POSIXct(paste(Sys.Date(),"00:00:00"))+600*(t-1))
        counter<-counter+1
        data[counter]<-mean(P1$P1[find&find2])
  }
}

data.matrix<-matrix(data,ncol=L)
plot(apply(data.matrix,2,mean,na.rm=TRUE))
image(data.matrix)

113 locations with 8038 observations

Simple spTimer model

library(spTimer)
GP<-spT.Gibbs(data~1,coords=cbind(locid$lon[-55],locid$lat[-55]))
newlat<-rep(seq(bbox[1],bbox[3],length=20),each=20)
newlon<-rep(seq(bbox[2],bbox[4],length=20),times=20)
newcoords<-data.frame("lat"=newlat,"lon"=newlon)
p.GP<-predict.spT(GP,newcoords=newcoords)


pred<-array(p.GP$Mean,c(20,144,20))

I<-5
newlat<-rep(seq(bbox[1],bbox[3],length=I),each=I)
newlon<-rep(seq(bbox[2],bbox[4],length=I),times=I)
newcoords<-cbind(newlat,newlon)
GPP<-spT.Gibbs(data~1,coords=100*cbind(locid$lon[-55],locid$lat[-55]), model="GPP", knots.coords=newcoords)
plot(GPP)

I<-12
newlat<-rep(seq(bbox[1],bbox[3],length=I)[-c(1,I)],each=I-2)
newlon<-rep(seq(bbox[2],bbox[4],length=I)[-c(1,I)],times=I-2)
newcoords<-100*cbind(newlat,newlon)
p.GPP<-predict.spT(GPP,newcoords=newcoords)

print(p.GPP)
pred<-array(p.GPP$Mean,c(144,20,20))
fields::image.plot(pred[1,,])
plot(apply(pred,1,mean))

Kriging (just spatial)

locations<-unique(locid)
L<-length(locations)
data<-rep(NA,L)
counter<-0
# bad coding ahead
for (j in 1:L)
{
    find2<-P1$locid==locations$id[j]
        counter<-counter+1
        data[counter]<-mean(P1$P1[find2])

}
data<-data.frame("P1"=data,"lat"=locid$lat,"lon"=locid$lon)

bbox2=array(bbox,c(2,2))
dimnames(bbox2)<-list(c("coords.x2","coords.x1"),c("min","max"))
locations2<-sp::SpatialPoints(cbind(data$lon,data$lat),CRS("+proj=longlat +datum=WGS84"),bbox=bbox2)
#coordinates(data)<-~lon+lat
P<-gstat(formula=data$P1~1,locations=locations2)

#v<-variogram(P)
#plot(v)
#fitP<-fit.variogram(v,vgm("Exp"))
#plot(v,fitP)
y<-rep(seq(bbox[2],bbox[4],length=100),each=100)
x<-rep(seq(bbox[1],bbox[3],length=100), times=100)
grid<-data.frame("lon"=x,"lat"=y,"NA"=NA)
grid<-sp::SpatialPoints(cbind(grid$lat,grid$lon),CRS("+proj=longlat +datum=WGS84"),bbox=bbox2)
z<-predict(P,grid)
pred<-t(array(z$var1.pred,c(100,100)))
fields::image.plot(seq(bbox[2],bbox[4],length=100),seq(bbox[1],bbox[3],length=100),pred, xlab="", ylab="")
points(locations2)

library(osmdata)
query <- opq(bbox='Munich, Germany')

query <- add_osm_feature(query, key='highway', value='secondary')
test<-osmdata_sp(query)
sp::plot(test$osm_lines, add=TRUE)

#P1.kriged <- krige(formula=P1 ~ 1, data=P, locations = ~lat+lon, P1, grid, model=fitP,bbox=bbox2)