code_03_test_Allan_data_v02.R

library(magrittr)
library(raster)

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following object is masked from 'package:magrittr':
## 
##     extract

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(FRK)

## 
## Attaching package: 'FRK'

## The following object is masked from 'package:raster':
## 
##     distance

library(spacetime)
library(sp)
library(rgdal)

## rgdal: version: 1.2-8, (SVN revision 663)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
##  Path to GDAL shared files: /usr/share/gdal/2.2
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.2-5

library(mapview)

## Loading required package: leaflet

library(fst)

# devtools::install_github("andrewzm/dggrids")

setwd("/home/michael/Dropbox/BGU/Itai_Kloog/p_46_FRK_Kriging_test/")

# Read data
dat = read.fst("01_raw/modX2009.fst")

# Rename
dat$lon = dat$Longitude
dat$lat = dat$Latitude
dat$Longitude = NULL
dat$Latitude = NULL

# Subset days
dat = dat[dat$day < (min(dat$day) + 11), ]

# Standard deviation of measurements (cannot be estimated from data in spatio-temporal kriging)
dat$std = 1

# Observations to 'spacetime' object
STObj = stConstruct(
  x = dat, # dataset
  space = c("lon","lat"), # spatial fields
  time = "day", # time field
  crs = CRS("+proj=longlat +ellps=sphere"), # CRS
  interval = TRUE # time reflects an interval
) 

###################################
# Plot
sp = as(STObj, "Spatial")
plot(sp)

###################################

# Prediction (BAU) grid
grid_BAUs = auto_BAUs(
  manifold = STsphere(), # space-time field on sphere
  data = STObj, 
  isea3h_res = 9,
  type = "hex", 
  tunit = "days" # time spacing in days
) 

## Loading required namespace: dggrids

## Loading required namespace: rgeos

grid_BAUs$fs = 1

###################################
# Plot
sp = as(grid_BAUs, "Spatial")
plot(sp)

###################################

# Basis functions
G_spatial = auto_basis(
  manifold = sphere(), # functions on sphere
  data = as(STObj,"Spatial"), # collapse time out
  nres = 3, # use three DGGRID resolutions
  prune = 15, # prune basis functions
  type = "bisquare", # bisquare basis functions
  subsamp = 2000, # use only 2000 data points for pruning
  isea3h_lo = 2
) 
G_temporal = local_basis(
  manifold = real_line(), # functions on real line
  loc = matrix(c(2,7,12)), # location parameter
  scale = rep(3,3), # scale parameter
  type = "Gaussian"
)
G_spacetime = TensorP(G_spatial,G_temporal)

# Model
f = error ~ 1 # formula for fixed effects
S = SRE(
  f = f, # formula
  data = list(STObj), # spatio-temporal object
  basis = G_spacetime, # space-time basis functions
  BAUs = grid_BAUs, # space-time BAUs
  est_error = FALSE, # do not estimate measurement error
  average_in_BAU = TRUE
) # average data that fall inside BAUs

## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.086 seconds
## Binned data in 0.087 seconds
## Binned data in 0.06 seconds
## Binned data in 0.078 seconds
## Binned data in 0.057 seconds
## Binned data in 0.054 seconds
## Binned data in 0.072 seconds
## Binned data in 0.058 seconds
## Binned data in 0.06 seconds
## Binned data in 0.076 seconds
## Binned data in 0.064 seconds

# Fit
S = SRE.fit(
  SRE_model = S, # SRE model
  n_EM = 2, # max. EM iterations
  tol = 0.01
) # convergence criteria

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |================================                                 |  50%
  |                                                                       
  |=================================================================| 100%
## Maximum EM iterations reached

# Predict
grid_BAUs = predict(S, obs_fs = TRUE) # predict only at select days

###################################
# Plot - Space-time
stplot(grid_BAUs[, , "mu"])

###################################
# Plot - Space only
grid_BAUs[, 3, "mu"] %>% spplot

###################################