# load required libraries
library(lgcp)
library(spatstat)
library(sp)
library(ncdf)
library(splancs)
#
# read data
ixyt<-read.table("../aegiss_ixyt.txt",header=T)
hampshire<-read.table("../aegiss_poly.txt",header=F)
#
par(pty="s")
polymap(hampshire)
pointmap(ixyt[,2:3],pch=19,cex=0.25,add=T)
# simplified boundary will speed up some computations
# (not recommended in real applications)
poly<-getpoly()
win<-owin(poly=poly)
# extract first six months’ data only, jitter to avoid duplicate
# locations and convert to kilometres
data<-ixyt[ixyt$t<=181,2:4]
data[,1]<-jitter(data[,1],amount=0.05)/1000
data[,2]<-jitter(data[,2],amount=0.05)/1000
tlim<-c(0,182)
xyt <-stppp(list(data=data,tlim=tlim,window=win))
par(mfrow=c(1,2),pty="s")
plot(xyt); plot(sort(xyt$t),1:length(xyt$t),pch=19,cex=0.5)
#
#
# estimate marginal spatial intensity (kernel estimator)
#
OW<-selectObsWindow(xyt,cellwidth=2) # grid-cells of width 2km
den <- lambdaEst(xyt,axes=TRUE)
plot(den)
#
# convert to spatial-at-risk object
#
sar <- spatialAtRisk(den)
#
# estimate temporal intensity (lowess of square-root-counts)
mut <- muEst(xyt,f=0.6)
plot(mut,ylim=c(0,10))
#
gin<-ginhomAverage(xyt,spatial.intensity=sar,
temporal.intensity=mut,rvals=0.025*(0:250))
plot(gin$r,gin$iso,type="l",xlab="r",ylab="g(r)")
kin <- KinhomAverage(xyt,spatial.intensity=sar,
temporal.intensity=mut,rvals=0.05*(0:250))
plot(kin$r,kin$iso-kin$theo,type="l",xlab="r",ylab="K(r)")
#
#
# minimum contrast parameter estimation: new code for next
# release of lgcp
#
source("minimumContrast.R")
minc <- minimum.contrast.spatiotemporal(xyt, model="exponential",
method="g",transform=log,spatial.dens=den,
temporal.intens=mut)