library(ggplot2)
library(bayesm)
library(zoo)
library(splines)
library(maptools)
library(plyr)
library(ggmap)
library(broom)
data("WaterQuality")
load("../../../Projects/2016 - SOB indicators/data/Maps/Maps.Rdata")
SOI = levels(as.factor(WQ_PPB$Species))[32]
names(WQ_PPB)
subWQ.df = subset(WQ_PPB, Species == SOI, select = c(Date,Measure,Site,Latitude, Longitude, Species))
head(subWQ.df)
range(subWQ.df$Date)
subWQ.df = subset(subWQ.df, Date < "2013-01-01", select = c(Date,Measure,Site,Latitude, Longitude, Species))
subWQ.df = subset(subWQ.df, Date > "2000-01-01", select = c(Date,Measure,Site,Latitude, Longitude, Species))
subWQ.df = subWQ.df[!is.na(subWQ.df$Measure),]
subWQ.df$Site <- tolower(sub("PORT PHILLIP BAY ","",subWQ.df$Site))
subWQ.df$Site = as.factor(subWQ.df$Site)
subWQ.df$Species = as.factor(subWQ.df$Species)
# plot #
bSP = bbox(SP)*matrix(c(0.999,1.0005,1.003,0.999),2,2, byrow=T)
rownames(bSP) <- c("lon", "lat")
bSP = asSpatialPolygonsBbox(bSP)
PPB_SHP <- spTransform(PPB_SHP, "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
Map_PPB <- gIntersection(PPB_SHP, bSP, byid = TRUE, drop_lower_td = TRUE)
SP = SpatialPointsDataFrame(coords = subWQ.df[,c("Longitude", "Latitude")], data = subWQ.df)
summaryPlotSPT(SP, bgmap = Map_PPB)
bg1 = Map_PPB
bg1.points = tidy(bg1)
#bg1@data$id = rownames(bg1@data)
bg2 = Habitat[Habitat@data$LC == "Subtidal Soft Sediment",]
bg1.df = bg1.points #join(bg1.points, bg1@data, by="id")
PPB_Map <- ggplot(bg1.df) + aes(long,lat, group = group) +
geom_polygon() +geom_path(color="white")
PPB_Map + theme(legend.position="bottom",legend.text=element_text(size=8),legend.key = element_blank() ) +
geom_point(aes(x = Longitude, y = Latitude, group =NULL, color =Site, size = Measure), data = SP@data) +
xlab("") + ylab("")
names(subWQ.df) <- c("date", "obs", "ID", "Latitude", "Longitude","Species")#c("ID" ,"date", "obs")
subWQ.df$ID <- as.factor(as.character(subWQ.df$ID))
subWQ.df$obs[subWQ.df$obs==0] <- NA
subWQ.df$obs <- log(subWQ.df$obs)
ggplot(data = subWQ.df, aes(x = date, y = obs, color = ID)) +geom_point() + facet_wrap(facets = ~ID) +
theme(legend.position = "bottom", legend.text=element_text(size=8),legend.key = element_blank() )
# Merging GMP and OMU as they seem to behave identically
subGroup.df = data.frame(Gr1 = rep(0.33, nlevels(subWQ.df$ID)),
Gr2 = rep(0.33, nlevels(subWQ.df$ID)),
Gr3 = rep(0.33, nlevels(subWQ.df$ID)), row.names = levels(subWQ.df$ID))
subCoords.df = unique(subWQ.df[,c("ID","Longitude", "Latitude")])
names(subCoords.df)[2:3] = tolower(names(subCoords.df)[2:3])
p.old= as.matrix(subGroup.df)
df = SPTMData(subWQ.df, agg = "quarter", seasonality = c("free"), debug=F,
method = "splines", degFree = 6)
p.old = p.old[match(levels(df$obs.data$ID),rownames(p.old)),]
y = df$obs.data
lu.long = data.frame(lu = as.factor(apply(subGroup.df[,1:ncol(subGroup.df)],1, which.max)), ID = as.character(rownames(subGroup.df)))
full = merge(y, lu.long, by="ID")
full = full[order(full$date),]
ggplot(data = full, aes(x = date, y = obs, col = lu)) + geom_point() + facet_wrap(~ID)
season = df$season
season.sites=season
season.lu = updateSeason(season.sites, p.old)
plot(season.lu, select = c(1:3))
theta.mean = rnorm(ncol(p.old), mean = rep(mean(y$obs, na.rm=T),ncol(p.old)),
sd = rep(mean(abs(y$obs), na.rm=T),ncol(p.old))/10)
theta.season = rnorm(ncol(season.lu), mean = 1, sd = 0.1)
theta.sigma = rnorm(ncol(p.old), mean = -0.5, sd = 0.05)
theta.0 = c(theta.mean,theta.season, theta.sigma)
y.p = meanFunc(theta.0, n.mixt=3, p.sites=p.old, season.sites = season.sites)
theta.mean = rnorm(ncol(p.old), mean = rep(mean(y$obs, na.rm=T),ncol(p.old)),
sd = rep(mean(abs(y$obs), na.rm=T),ncol(p.old))/10)
theta.season = rnorm(ncol(df$basis)*ncol(p.old), mean = 1, sd = 0.5)
theta.sigma = rnorm(ncol(p.old), mean = -0.5, sd = 0.05)
theta.0 = c(theta.mean,theta.season, theta.sigma)
y.p = meanFuncSplines(theta.0, n.mixt=3, basis = df$basis)
SiteID = unique(df$obs.data$ID)
par(mfrow=c(3,3), mar=c(2,2,2,1))
for(i in 1:length(SiteID)){
obs = y$obs[y$ID %in% SiteID[i]]
idx.k = which.max(p.old[i,])
pred= y.p[,idx.k]
dates = y$date[y$ID %in% SiteID[i]]
ci = cbind(pred - 2*exp(theta.0[(length(theta.0) - ncol(p.old) + 1):length(theta.0)][idx.k]),pred + 2*exp(theta.0[(length(theta.0) - ncol(p.old) + 1):length(theta.0)][idx.k]) )
plot(dates,obs,main = SiteID[i],
type="l", ylim = range(ci, na.rm=T))
polygon(x=c(rev(dates), (dates)), y = c(rev(ci[,1]), (ci[,2])),
border=NA, col = "grey", density = 100)
points(dates,pred, type="l", col="red")
points(y$date[y$ID %in% SiteID[i]],y$obs[y$ID %in% SiteID[i]])
}
subCoords.df = subCoords.df[match(levels(df$obs.data$ID),subCoords.df$ID),]
proj.xy = project(cbind(subCoords.df$longitude,subCoords.df$latitude), "+proj=utm +zone=51 ellps=WGS84")/1000
subCoords.df$x = proj.xy[,1]
subCoords.df$y = proj.xy[,2]
sptmod = SPTModel(df, p.old, proj.xy)
resEM = estimEM(sptmod, theta.init = NULL, print.res = T, debug=F)
theta.EM = resEM$Result$theta
p.EM = resEM$Result$latent
rownames(p.EM) = rownames(p.old)
df.p.EM = as.data.frame(p.EM)
df.p.EM$Site = rownames(p.EM)
df.p.EM$Group = as.factor(apply(df.p.EM[,1:3],1,which.max))
plot(resEM, type = "residuals")
plot(resEM, type = "tempBasis")
plot(resEM, type = "prediction")
colnames(p.EM)<-colnames(p.old)
season.lu = updateSeason(season.sites, p.EM)
season.lu[is.na(season.lu)] <- NA
par(mar = c(5,4,4,1))
plot(season.lu, select = 1:3)
SP@data = merge(SP@data, df.p.EM, by = "Site")
PPB_Map + theme(legend.position="bottom",legend.text=element_text(size=8),legend.key = element_blank() ) +
geom_point(aes(x = Longitude, y = Latitude, group =NULL, color =Group, size = Measure), data = SP@data) +
xlab("") + ylab("")
theta.EM[(length(theta.EM) - ncol(p.EM) + 1):length(theta.EM)] = exp(theta.EM[(length(theta.EM) - ncol(p.EM) + 1):length(theta.EM)])
meanFunc(theta.EM, n.mixt = 5, p.EM, season.sites)
resGibbs = estimGibbs(sptmod, method = 'simple', theta.init = theta.EM, print.res = T,
N.run = 1000, debug=FALSE)
plot(resGibbs)
par(mar=c(0,0,0,0))
diagPlot(resGibbs, par = "theta",
main = "", breaks = 50)
diagPlot(resGibbs, par = "paired-theta",
main = "", breaks = 50)
diagPlot(resGibbs, par = "tempBasis")
theta.Gibbs = resGibbs$Result$theta
p.Gibbs = resGibbs$Result$latent
season.lu = updateSeason(season.sites, p.Gibbs)
season.lu[is.na(season.lu)] <- NA
plot(season.lu, select = 1:5)
SpPt = SpatialPointsDataFrame(coords = as.matrix(subCoords.df[,2:3]), data = data.frame(Comp = apply(p.Gibbs, 1, which.max)), proj4string = CRS("+proj=longlat +ellps=GRS80 +no_defs"))
library(leaflet)
m <- leaflet() %>%
addTiles() %>%
addPolygons(data = CatchmentMap[which(CatchmentMap@data$Site_Code %in% SiteID),],
stroke = FALSE, fillOpacity = 0.2, smoothFactor = 0.5) %>%
addCircles(data = SpPt, color = c("green", "red", "blue","orange","pink")[SpPt@data$Comp], radius=rep(400,nrow(subCoords.df)))
# addMarkers(data = SpPt, popup = as.character(unlist(Coords[,1])))
m
resGibbsSpat = estimGibbs(sptmod, method = 'spat', theta.init = theta.EM, print.res = T,
N.run = 15000, beta.range = c(0,3))
plot(resGibbsSpat)
pred = predict(resGibbs)
diagPlot(resGibbsSpat, par = "theta",
main = "", breaks = 50)
diagPlot(resGibbsSpat, par = "paired-theta",
main = "", breaks = 50)
diagPlot(resGibbsSpat, par = "beta",
main = "", breaks = 50)
p.GibbsSpat = resGibbsSpat$Result$latent
season.lu = updateSeason(season.sites, p.GibbsSpat)
plot(season.lu)
SpPt = SpatialPointsDataFrame(coords = as.matrix(subCoords.df[,2:3]), data = data.frame(Comp = apply(p.GibbsSpat, 1, which.max)), proj4string = CRS("+proj=longlat +ellps=GRS80 +no_defs"))
library(leaflet)
m <- leaflet() %>%
addTiles() %>%
addPolygons(data = CatchmentMap[which(CatchmentMap@data$Site_Code %in% SiteID),],
stroke = FALSE, fillOpacity = 0.2, smoothFactor = 0.5) %>%
addCircles(data = SpPt, color = c("green", "red", "blue","orange")[SpPt@data$Comp], radius=rep(200,nrow(subCoords.df)))
# addMarkers(data = SpPt, popup = as.character(unlist(Coords[,1])))
m
par(mfrow=c(1,2), mar = c(1,1,1,1))
plot(CatchmentMap[which(CatchmentMap@data$Site_Code %in% SiteID),])
points(subCoords.df$longitude, subCoords.df$latitude)
compCircle(subCoords.df$longitude, subCoords.df$latitude, p.old, color = rainbow(6), radius = 0.015)
text(subCoords.df$longitude, subCoords.df$latitude, subCoords.df$ID, cex=.7)
plot(CatchmentMap[which(CatchmentMap@data$Site_Code %in% SiteID),])
points(subCoords.df$longitude, subCoords.df$latitude)
compCircle(subCoords.df$longitude, subCoords.df$latitude, p.Gibbs, color = rainbow(6), radius = 0.015)
text(subCoords.df$longitude, subCoords.df$latitude, subCoords.df$ID, cex=.7)
# Voronoi Cells
cov = subCoords.df
library(deldir)
library(rgeos)
par(mfrow = c(1,3), mar = c(5,4,4,1))
plot(cov$longitude, cov$latitude, xlab = "Longitude", ylab = "Latitude")
dd = deldir(cov$longitude, cov$latitude)
tile.dd = tile.list(dd)
col.fill = rainbow(8, alpha = 0.2)[apply(p.old, 1, which.max)]
plot(tile.dd, add=T, fillcol = col.fill, border = NA, close = FALSE)
plot(CatchmentMap[which(CatchmentMap@data$Site_Code %in% SiteID),], add=T)
plot(CatchmentMap, add=T)
plot(cov$longitude, cov$latitude)
col.fill = rainbow(8, alpha = 0.2)[apply(p.Gibbs, 1, which.max)]
plot(tile.dd, add=T, fillcol = col.fill, border = NA, close = FALSE)
plot(CatchmentMap, add=T)
plot(cov$longitude, cov$latitude)
col.fill = rainbow(8, alpha = 0.2)[apply(p.GibbsSpat, 1, which.max)]
plot(tile.dd, add=T, fillcol = col.fill, border = NA, close = FALSE)
plot(CatchmentMap, add=T)
library(spatstat)
library(maptools)
df = data.frame(x = cov$longitude, y = cov$latitude)
W <- ripras(df, shape="rectangle")
W <- owin(range(cov$longitude)*c(0.999,1.001), range(cov$latitude)*c(1.001,0.999))
# create point pattern
X <- as.ppp(df, W=W)
plot(X)
# compute Dirichlet triangles
Y <- dirichlet(X)
plot(Y, add=T)
# convert
Z <- as(Y, "SpatialPolygons")
Xp <- as(X, "SpatialPoints")
par(mfrow= c(1,1), mar = c(5,4,4,1))
plot(Z)
plot(subCatchmentMap, add=T)
CatMap = CatchmentMap[which(colSums(gIntersects(CatchmentMap, Xp, byid = TRUE)) > 0),]
CatMap2 = CatchmentMap[which(CatchmentMap@data$CATNAME == "Onkaparinga River" & CatchmentMap@data$Site_Code %in% SiteID),]
plot(CatMap)
ssCM = gIntersection(CatMap,Z, byid = TRUE, drop_lower_td=TRUE)
IDs = NULL
for(i in 1:length(ssCM@polygons)){
IDs = rbind(IDs, c(unlist(strsplit(ssCM@polygons[[i]]@ID, " "))))
}
palette <- colorRampPalette(c("green","red"))(7)
palette[sort(resGibbsSpat$Result$theta[1:7], index.return=T)$ix] = palette
col.fill= palette[apply(p.GibbsSpat, 1, which.max)]
par(mfrow=c(1,1))
plot(CatMap2, ylim = c(-35.20,-34.87))
axis(1);axis(2)
plot(ssCM, col = col.fill[as.numeric(IDs[,2])], border=NA, add=T)
plot(CatMap, add=T)
plot(X, add=T, cex = .5)
legend(x = 138.43,y=-34.90,legend = names(subLandUse.df)[sort(resGibbsSpat$Result$theta[1:7], index.return=T)$ix], fill = palette[sort(resGibbsSpat$Result$theta[1:7], index.return=T)$ix])
pdf(file = "../MapNitrateLevels.pdf", width = 10, height = 8)
par(mfrow=c(1,1))
plot(CatchmentMap, xlim = c(138.4,139.0), ylim = c(-35.2,-34.93))
plot(CatMap, add=T)
axis(1);axis(2)
plot(ssCM, col = col.fill[as.numeric(IDs[,2])], border=NA, add=T)
plot(CatMap, add=T)
plot(X, add=T, cex = .5)
legend(x = 138.41,y=-34.96,cex = 0.85,
legend = names(subLandUse.df)[sort(resGibbs$Result$theta[1:7], index.return=T)$ix],
fill = palette[sort(resGibbs$Result$theta[1:7], index.return=T)$ix])
dev.off()
require(graphics)
library(splines)
library(MASS)
splineDesign(knots = 1:10, x = 4:7)
splineDesign(knots = 1:10, x = 4:7, deriv = 1)
## visualize band structure
Matrix::drop0(zapsmall(6*splineDesign(knots = 1:40, x = 4:37, sparse = TRUE)))
knots <- seq(1,10,.85) # 10 => 11-4 = 7 Basis splines
x <- seq(min(knots)+1.5, max(knots)-1.5, length.out = 501)
bb <- splineDesign(knots, x = x, outer.ok = TRUE)
dataT = df$obs.data
y=dataT$obs[!is.na(dataT$obs)] - mean(dataT$obs[!is.na(dataT$obs)])
x=dataT$date[!is.na(dataT$obs)]
plot(x,y, type="p")
bbW = bs(as.numeric(format(x, "%W")),df = 51)
bbS = bs(as.numeric(format(x, "%m")),df = 11)
bbT = bs(as.numeric(format(x, "%Y")),df = diff(range(as.numeric(format(x, "%Y")))))
bbF = cbind(bbW,bbS,bbT)
w = ginv(t(bbF) %*% bbF) %*% t(bbF) %*% matrix(y, ncol=1)
var.w = ginv(t(bbF) %*% bbF) * mean((y-bbF %*% w)^2)
plot(x,(y), type="p")
points(x, bbF %*% w, type= "l", lwd=2, col="red")
points(x, bbF %*% w + 1.96*sqrt(diag((bbF) %*% var.w %*% t(bbF))), type="l", lwd = 0.5, col="red")
points(x, bbF %*% w - 1.96*sqrt(diag((bbF) %*% var.w %*% t(bbF))), type="l", lwd = 0.5, col="red")
idx.w = 1:ncol(bbW)
df.week = unique(cbind(as.numeric(format(x, "%W")),bbW %*% w[idx.w],
bbW %*% w[idx.w] - 1.96 * sqrt(diag((bbW) %*% var.w[idx.w,idx.w] %*% t(bbW))),
bbW %*% w[idx.w] + 1.96 * sqrt(diag((bbW) %*% var.w[idx.w,idx.w] %*% t(bbW)))))
plot(df.week[,1],df.week[,2], type="p", ylim = range(df.week[,2:4]))
arrows(df.week[,1], df.week[,3], df.week[,1], df.week[,4], code=0)
abline(h=0)
idx.m = (ncol(bbW)+1):(ncol(bbW)+ncol(bbS))
df.month = unique(cbind(as.numeric(format(x, "%m")),bbS %*% w[idx.m],
bbS %*% w[idx.m] - 1.96 * sqrt(diag((bbS) %*% var.w[idx.m,idx.m] %*% t(bbS))),
bbS %*% w[idx.m] + 1.96 * sqrt(diag((bbS) %*% var.w[idx.m,idx.m] %*% t(bbS)))))
plot(df.month[,1],df.month[,2], type="l", ylim = range(df.month[,2:4]))
arrows(df.month[,1], df.month[,3], df.month[,1], df.month[,4], code=0)
abline(h=0)
idx.y = ((ncol(bbW)+ncol(bbS))+1):length(w)
df.year = unique(cbind(as.numeric(format(x, "%Y")),bbT %*% w[idx.y],
bbT %*% w[idx.y] - 1.96 * sqrt(diag((bbT) %*% var.w[idx.y,idx.y] %*% t(bbT))),
bbT %*% w[idx.y] + 1.96 * sqrt(diag((bbT) %*% var.w[idx.y,idx.y] %*% t(bbT)))))
plot(df.year[,1],df.year[,2], type="l", ylim = range(df.year[,2:4]))
arrows(df.year[,1], df.year[,3], df.year[,1], df.year[,4], code=0)
abline(h=0)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.