SSN
Spatial Modeling on Stream Networks

Package index
Search the SSN package
Vignettes
Functions
166
Source code
80
Man pages
54
Home
/
CRAN
/
SSN
/
inst/doc/SSN.R

inst/doc/SSN.R
In SSN: Spatial Modeling on Stream Networks 

## ----results='hide',echo=FALSE, cache=FALSE------------------------------
## Code to run in the background
set.seed(210112)

## ------------------------------------------------------------------------
library("SSN")

## ------------------------------------------------------------------------
file.copy(system.file("lsndata/MiddleFork04.ssn", package = "SSN"),
  to = tempdir(), recursive = TRUE, copy.mode = FALSE)
setwd(tempdir())

## ----eval=FALSE----------------------------------------------------------
## importSSN(Path, predpts = NULL, o.write = FALSE)

## ------------------------------------------------------------------------
mf04p <- importSSN("./MiddleFork04.ssn",
   predpts = "pred1km")

## ------------------------------------------------------------------------
mf04p <- importPredpts(mf04p, "Knapp", "ssn")
mf04p <- importPredpts(mf04p, "CapeHorn", "ssn")

## ----eval=FALSE----------------------------------------------------------
## additive.function(mf04p, VarName, afvName)

## ------------------------------------------------------------------------
names(mf04p@data)

## ------------------------------------------------------------------------
head(mf04p@data[, c("h2oAreaKm2", "afvArea")])

## ------------------------------------------------------------------------
mf04p <- additive.function(mf04p, "h2oAreaKm2", "computed.afv")

## ------------------------------------------------------------------------
names(mf04p@data)

## ------------------------------------------------------------------------
head(mf04p@data[, c("h2oAreaKm2",
   "afvArea", "computed.afv")])
head(getSSNdata.frame(mf04p)[, c("afvArea", "computed.afv")])

## ------------------------------------------------------------------------
createDistMat(mf04p, predpts = "Knapp", o.write = TRUE,
	amongpreds = TRUE)
createDistMat(mf04p, predpts = "CapeHorn", o.write = TRUE,
	amongpreds = TRUE)

## ------------------------------------------------------------------------
distObs <- getStreamDistMat(mf04p)
str(distObs)
distObs$dist.net1[1:5,1:5]

## ------------------------------------------------------------------------
strDistNet2 <- distObs$dist.net2 + t(distObs$dist.net2)
strDistNet2[5:10,5:10]

## ------------------------------------------------------------------------
distPred1km <- getStreamDistMat(mf04p, Name = "pred1km")
str(distPred1km)
distPred1km$dist.net1.a[1:5,1:5]

## ------------------------------------------------------------------------
createDistMat(mf04p, predpts = "CapeHorn", o.write = TRUE, 
  amongpreds = TRUE)
distCape <- getStreamDistMat(mf04p, Name = "CapeHorn")
str(distCape)
distCape$dist.net2[1:5,1:5]

## ------------------------------------------------------------------------
names(mf04p)

## ----LoadData,include=FALSE,results='hide',cache=FALSE-------------------
plot(mf04p, lwdLineCol = "afvArea", lwdLineEx = 10, lineCol = "blue",
   pch = 19, xlab = "x-coordinate (m)", ylab = "y-coordinate (m)",
   asp = 1)

## ----eval=FALSE----------------------------------------------------------
## plot(mf04p, lwdLineCol = "afvArea", lwdLineEx = 10, lineCol = "blue",
##    pch = 19, xlab = "x-coordinate (m)", ylab = "y-coordinate (m)",
##    asp = 1)

## ----plotSpatialStreamNetwork,include=FALSE,results='hide', cache=FALSE----
brks <- plot(mf04p, "Summer_mn", lwdLineCol = "afvArea",
   lwdLineEx = 15, lineCol = "black", xlab =  "x-coordinate" ,
   ylab =  "y-coordinate", asp=1 )

## ----eval=FALSE----------------------------------------------------------
## brks <- plot(mf04p, "Summer_mn", lwdLineCol = "afvArea",
##    lwdLineEx = 15, lineCol = "black", xlab =  "x-coordinate" ,
##    ylab =  "y-coordinate", asp=1 )

## ----plotUsingSP, include=FALSE, results='hide', cache=FALSE-------------
#plot the stream lines
plot(as.SpatialLines(mf04p), col = "blue")
# add the observed locations with size proportional 
# to mean summer temperature
plot(as.SpatialPoints(mf04p), pch = 19, 
  cex = as.SpatialPointsDataFrame(mf04p)$Summer_mn/9 , add = TRUE)
# add the prediction locations on the 1 km spacing
plot(as.SpatialPoints(mf04p, data = "pred1km"), cex = 1.5, add = TRUE)
# add the dense set of points for block prediction on Knapp segment
plot(as.SpatialPoints(mf04p, data = "Knapp"), pch = 19, cex = 0.3, 
  col = "red", add = TRUE)

## ----eval=FALSE----------------------------------------------------------
## #plot the stream lines
## plot(as.SpatialLines(mf04p), col = "blue")
## # add the observed locations with size proportional
## # to mean summer temperature
## plot(as.SpatialPoints(mf04p), pch = 19,
##   cex = as.SpatialPointsDataFrame(mf04p)$Summer_mn/9 , add = TRUE)
## # add the prediction locations on the 1 km spacing
## plot(as.SpatialPoints(mf04p, data = "pred1km"), cex = 1.5, add = TRUE)
## # add the dense set of points for block prediction on Knapp segment
## plot(as.SpatialPoints(mf04p, data = "Knapp"), pch = 19, cex = 0.3,
##   col = "red", add = TRUE)

## ----Torgegram, include=FALSE,cache=FALSE--------------------------------
mf04.Torg <- Torgegram(mf04p, "Summer_mn", nlag = 20, maxlag = 50000)
plot(mf04.Torg)

## ----eval=FALSE----------------------------------------------------------
## mf04.Torg <- Torgegram(mf04p, "Summer_mn", nlag = 20, maxlag = 50000)
## plot(mf04.Torg)

## ----GauModel0-----------------------------------------------------------
mf04.glmssn0 <- glmssn(Summer_mn ~ ELEV_DEM + SLOPE, mf04p,
   CorModels = NULL, use.nugget = TRUE)
summary(mf04.glmssn0)

## ----lm0-----------------------------------------------------------------
summary(lm(Summer_mn ~ ELEV_DEM + SLOPE, getSSNdata.frame(mf04p)))

## ----GauModel1, cache = TRUE---------------------------------------------
mf04.glmssn1 <- glmssn(Summer_mn ~ ELEV_DEM + SLOPE, mf04p,
   CorModels = c("Exponential.tailup", "Exponential.taildown",
      "Exponential.Euclid"), addfunccol = "afvArea")
summary(mf04.glmssn1)

## ----BinModel1, cache = TRUE---------------------------------------------
mf04.glmssnBin <- glmssn(MaxOver20 ~ ELEV_DEM + SLOPE, mf04p,
  CorModels = c("Mariah.tailup", "Spherical.taildown"),
  family = "binomial", addfunccol = "afvArea")
summary(mf04.glmssnBin)

## ----PoiModel1, cache = TRUE---------------------------------------------
mf04.glmssnPoi <- glmssn(C16 ~ ELEV_DEM + SLOPE, mf04p,
  CorModels = c("LinearSill.tailup", "LinearSill.taildown"),
  family = "poisson", addfunccol = "afvArea")
summary(mf04.glmssnPoi)

## ----Model1, include=FALSE, cache=FALSE----------------------------------
mf04.resid1 <- residuals(mf04.glmssn1)
names( getSSNdata.frame(mf04.resid1) )
plot(mf04.resid1)

## ------------------------------------------------------------------------
mf04.resid1 <- residuals(mf04.glmssn1)
names( getSSNdata.frame(mf04.resid1) )

## ----eval=FALSE----------------------------------------------------------
## plot(mf04.resid1)

## ----ResidHist,include=FALSE,cache=FALSE---------------------------------
par(mfrow = c(1, 2))
hist(mf04.resid1)
hist(mf04p, "Summer_mn")

## ----eval=FALSE----------------------------------------------------------
## par(mfrow = c(1, 2))
## hist(mf04.resid1)
## hist(mf04p, "Summer_mn")

## ------------------------------------------------------------------------
ObsDFr <- getSSNdata.frame(mf04.resid1)
ObsDF <- getSSNdata.frame(mf04p)
indOutlier <- ObsDFr["_resid_"] < -3
ObsDF[indOutlier, "Summer_mn"] <- NA
mf04c <- putSSNdata.frame(ObsDF, mf04p)

## ----cache=FALSE---------------------------------------------------------
mf04c.glmssn0 <- glmssn(Summer_mn ~ ELEV_DEM + SLOPE, mf04c,
   CorModels = c("Exponential.tailup", "Exponential.taildown",
   "Exponential.Euclid"), addfunccol = "afvArea", EstMeth = "ML")
summary(mf04c.glmssn0)

## ----cache=FALSE---------------------------------------------------------
mf04c.glmssn1 <- glmssn(Summer_mn ~ ELEV_DEM, mf04c,
   CorModels = c("Exponential.tailup", "Exponential.taildown"),
   addfunccol = "afvArea", EstMeth = "ML")
summary(mf04c.glmssn1)

## ----LOOCV, include=FALSE------------------------------------------------
cv.out <- CrossValidationSSN(mf04c.glmssn1)
par(mfrow = c(1, 2))
plot(mf04c.glmssn1$sampinfo$z,
   cv.out[, "cv.pred"], pch = 19,
   xlab = "Observed Data", ylab = "LOOCV Prediction")
abline(0, 1)
plot( na.omit( getSSNdata.frame(mf04c)[, "Summer_mn"]),
   cv.out[, "cv.se"], pch = 19,
   xlab = "Observed Data", ylab = "LOOCV Prediction SE")

## ----eval=FALSE----------------------------------------------------------
## cv.out <- CrossValidationSSN(mf04c.glmssn1)
## par(mfrow = c(1, 2))
## plot(mf04c.glmssn1$sampinfo$z,
##    cv.out[, "cv.pred"], pch = 19,
##    xlab = "Observed Data", ylab = "LOOCV Prediction")
## abline(0, 1)
## plot( na.omit( getSSNdata.frame(mf04c)[, "Summer_mn"]),
##    cv.out[, "cv.se"], pch = 19,
##    xlab = "Observed Data", ylab = "LOOCV Prediction SE")

## ----LOOCVSummary--------------------------------------------------------
CrossValidationStatsSSN(mf04c.glmssn1)

## ------------------------------------------------------------------------
GR2(mf04c.glmssn1)
varcomp(mf04c.glmssn1)

## ------------------------------------------------------------------------
AIC(mf04c.glmssn0)
AIC(mf04c.glmssn1)

## ------------------------------------------------------------------------
mf04c.glmssn1 <- glmssn(Summer_mn ~ ELEV_DEM, mf04c,
   CorModels = c("Exponential.tailup", "Exponential.taildown"),
   addfunccol = "afvArea")
mf04c.glmssn2 <- glmssn(Summer_mn ~ ELEV_DEM,  mf04c,
   CorModels = c("LinearSill.tailup", "Mariah.taildown"),
   addfunccol = "afvArea")
mf04c.glmssn3 <- glmssn(Summer_mn ~ ELEV_DEM , mf04c,
   CorModels =  c("Mariah.tailup", "LinearSill.taildown"),
   addfunccol = "afvArea")
mf04c.glmssn4 <- glmssn(Summer_mn ~ ELEV_DEM, mf04c,
   CorModels = c("Spherical.tailup", "Spherical.taildown"),
   addfunccol = "afvArea")
mf04c.glmssn5 <- glmssn(Summer_mn ~ ELEV_DEM, mf04c,
   CorModels = "Exponential.Euclid",
   addfunccol = "afvArea")

## ------------------------------------------------------------------------
options(digits = 4)
InfoCritCompare(list(mf04c.glmssn1, mf04c.glmssn2,
   mf04c.glmssn3, mf04c.glmssn4, mf04c.glmssn5))
options(digits = 7)

## ------------------------------------------------------------------------
summary(mf04c.glmssn2)

## ----Residuals,include=FALSE,cache=FALSE---------------------------------
mf04c.resid2 <- residuals(mf04c.glmssn2,
    cross.validation = TRUE)
mf04c.resid2.cv.std <-
    getSSNdata.frame(mf04c.resid2)[, "_resid.crossv_"] /
    getSSNdata.frame(mf04c.resid2)[, "_CrossValStdErr_"]
hist(mf04c.resid2.cv.std)

## ----eval=FALSE----------------------------------------------------------
## mf04c.resid2 <- residuals(mf04c.glmssn2,
##     cross.validation = TRUE)
## mf04c.resid2.cv.std <-
##     getSSNdata.frame(mf04c.resid2)[, "_resid.crossv_"] /
##     getSSNdata.frame(mf04c.resid2)[, "_CrossValStdErr_"]
## hist(mf04c.resid2.cv.std)

## ----TorgRes,include=FALSE,cache=FALSE-----------------------------------
plot(Torgegram(mf04c.resid2, "_resid_", nlag = 8, maxlag = 25000))

## ----eval=FALSE----------------------------------------------------------
## plot(Torgegram(mf04c.resid2, "_resid_", nlag = 8, maxlag = 25000))

## ----Preds1,include=FALSE,cache=FALSE------------------------------------
mf04c.pred1km <- predict(mf04c.glmssn4, "pred1km")
plot(mf04c.pred1km, SEcex.max = 1, SEcex.min = .5/3*2,
     breaktype = "user", brks = brks)

## ----eval=FALSE----------------------------------------------------------
## mf04c.pred1km <- predict(mf04c.glmssn4, "pred1km")
## plot(mf04c.pred1km, SEcex.max = 1, SEcex.min = .5/3*2,
##      breaktype = "user", brks = brks)

## ----Preds2,include=FALSE------------------------------------------------
plot(mf04c, "Summer_mn", pch = 1, cex = 3,
   xlab = "x-coordinate", ylab = "y-coordinate",
   xlim = c(-1511000,-1500000), ylim = c(2525000,2535000))
mf04c.glmssn4.Knapp <- predict(mf04c.glmssn4, "Knapp")
plot(mf04c.glmssn4.Knapp, "Summer_mn", add = TRUE,
   xlim = c(-1511000,-1500000), ylim = c(2525000,2535000))

## ----eval=FALSE----------------------------------------------------------
## plot(mf04c, "Summer_mn", pch = 1, cex = 3,
##    xlab = "x-coordinate", ylab = "y-coordinate",
##    xlim = c(-1511000,-1500000), ylim = c(2525000,2535000))
## mf04c.glmssn4.Knapp <- predict(mf04c.glmssn4, "Knapp")
## plot(mf04c.glmssn4.Knapp, "Summer_mn", add = TRUE,
##    xlim = c(-1511000,-1500000), ylim = c(2525000,2535000))

## ------------------------------------------------------------------------
mf04c.glmssn4.BPKnapp <- BlockPredict(mf04c.glmssn4, "Knapp")
mf04c.glmssn4.BPKnapp

## ------------------------------------------------------------------------
mf04c.glmssn4.BPCapeHorn <- BlockPredict(mf04c.glmssn4, "CapeHorn")
mf04c.glmssn4.BPCapeHorn

## ------------------------------------------------------------------------
mf04c.missingobs <- predict(mf04c.glmssn4, "_MissingObs_")
getPreds(mf04c.missingobs, pred.type = "pred")
with(getSSNdata.frame(mf04p), Summer_mn[pid==29])

## ----eval=FALSE----------------------------------------------------------
## createSSN(n, obsDesign, predDesign = noPoints, path,
##    importToR = FALSE, treeFunction = igraphKamadaKawai)

## ----SimIterative, include=FALSE, cache=FALSE----------------------------
set.seed(12)
iterative.ssn <- createSSN(n = c(30, 10),
   obsDesign = binomialDesign(c(10,10)),
   importToR = TRUE, path = "./SimIterative.ssn",
   treeFunction = iterativeTreeLayout)
plot(iterative.ssn, lwdLineCol = "addfunccol", lwdLineEx = 8,
   lineCol = "blue", cex = 2, xlab = "x-coordinate",
   ylab = "y-coordinate", pch = 1)

## ----eval=FALSE, cache=FALSE---------------------------------------------
## set.seed(12)
## iterative.ssn <- createSSN(n = c(30, 10),
##    obsDesign = binomialDesign(c(10,10)),
##    importToR = TRUE, path = "./SimIterative.ssn",
##    treeFunction = iterativeTreeLayout)
## plot(iterative.ssn, lwdLineCol = "addfunccol", lwdLineEx = 8,
##    lineCol = "blue", cex = 2, xlab = "x-coordinate",
##    ylab = "y-coordinate", pch = 1)

## ----SimSSN1, include=FALSE, cache=FALSE---------------------------------
set.seed(101)
raw.ssn <- createSSN(n = c(10, 10, 10),
   obsDesign = binomialDesign(c(40, 40, 40)),
   predDesign = systematicDesign(c(0.2, 0.4, 0.8)), importToR = TRUE,
   path = "./raw.ssn")
plot(raw.ssn, lwdLineCol = "addfunccol", lwdLineEx = 8,
   lineCol = "blue", cex = 2, xlab = "x-coordinate",
   ylab = "y-coordinate", pch = 1)
plot(raw.ssn, PredPointsID = "preds", add = TRUE, cex = .5, pch = 19,
   col = "green")

## ----eval=FALSE, cache=FALSE---------------------------------------------
## set.seed(101)
## raw.ssn <- createSSN(n = c(10, 10, 10),
##    obsDesign = binomialDesign(c(40, 40, 40)),
##    predDesign = systematicDesign(c(0.2, 0.4, 0.8)), importToR = TRUE,
##    path = "./raw.ssn")
## plot(raw.ssn, lwdLineCol = "addfunccol", lwdLineEx = 8,
##    lineCol = "blue", cex = 2, xlab = "x-coordinate",
##    ylab = "y-coordinate", pch = 1)
## plot(raw.ssn, PredPointsID = "preds", add = TRUE, cex = .5, pch = 19,
##    col = "green")

## ----SimHardcore, include=FALSE, cache=FALSE-----------------------------
set.seed(13)
hardcore.ssn <- createSSN(n = c(10, 10),
   obsDesign = hardCoreDesign(c(200, 200), c(0.2, 0.4)),
   importToR = TRUE, path = "./SimHardcore.ssn")
plot(hardcore.ssn, lwdLineCol = "addfunccol", lwdLineEx = 8,
   lineCol = "blue", cex = 2, xlab = "x-coordinate",
   ylab = "y-coordinate", pch = 1)
plot(hardcore.ssn, PredPointsID = NULL, add = TRUE, cex = .5,
   pch = 19, col = "green")

## ----eval=FALSE----------------------------------------------------------
## set.seed(13)
## hardcore.ssn <- createSSN(n = c(10, 10),
##    obsDesign = hardCoreDesign(c(200, 200), c(0.2, 0.4)),
##    importToR = TRUE, path = "./SimHardcore.ssn")
## plot(hardcore.ssn, lwdLineCol = "addfunccol", lwdLineEx = 8,
##    lineCol = "blue", cex = 2, xlab = "x-coordinate",
##    ylab = "y-coordinate", pch = 1)
## plot(hardcore.ssn, PredPointsID = NULL, add = TRUE, cex = .5,
##    pch = 19, col = "green")

## ----cache=FALSE---------------------------------------------------------
createDistMat(raw.ssn, "preds", o.write=TRUE, amongpred = TRUE)

## ------------------------------------------------------------------------
rawDFobs <- getSSNdata.frame(raw.ssn, "Obs")
rawDFpred <- getSSNdata.frame(raw.ssn, "preds")

## ------------------------------------------------------------------------
rawDFobs[,"X1"] <- rnorm(length(rawDFobs[,1]))
rawDFpred[,"X1"] <- rnorm(length(rawDFpred[,1]))
rawDFobs[,"X2"] <- rnorm(length(rawDFobs[,1]))
rawDFpred[,"X2"] <- rnorm(length(rawDFpred[,1]))

## ------------------------------------------------------------------------
rawDFobs[,"F1"] <- as.factor(sample.int(4,length(rawDFobs[,1]),
   replace = TRUE))
rawDFpred[,"F1"] <- as.factor(sample.int(4,length(rawDFpred[,1]),
   replace = TRUE))

## ------------------------------------------------------------------------
rawDFobs[,"RE1"] <- as.factor(sample(1:3,length(rawDFobs[,1]),
   replace = TRUE))
rawDFobs[,"RE2"] <- as.factor(sample(1:4,length(rawDFobs[,1]),
   replace = TRUE))
rawDFpred[,"RE1"] <- as.factor(sample(1:3,length(rawDFpred[,1]),
   replace = TRUE))
rawDFpred[,"RE2"] <- as.factor(sample(1:4,length(rawDFpred[,1]),
   replace = TRUE))

## ------------------------------------------------------------------------
names(rawDFobs)
names(rawDFpred)

## ----cache=FALSE---------------------------------------------------------
set.seed(102)
sim.out <- SimulateOnSSN(raw.ssn, ObsSimDF = rawDFobs,
   PredSimDF = rawDFpred, PredID = "preds",
   formula = ~ X1 + X2 + F1, coefficients = c(10,1,0,-2,0,2),
   CorModels = c("LinearSill.tailup", "Mariah.taildown",
   "Exponential.Euclid", "RE1", "RE2"), use.nugget = TRUE,
   CorParms = c(3, 10, 2, 10, 1, 5, 1, .5, .1),
   addfunccol = "addfunccol")

## ------------------------------------------------------------------------
with(rawDFobs, colnames(model.matrix( ~ X1 + X2 + F1)))

## ------------------------------------------------------------------------
sim.out$FixedEffects
sim.out$CorParms

## ------------------------------------------------------------------------
sim.ssn <- sim.out$ssn.object

## ----SimSSN2,include=FALSE,results='hide',cache=FALSE--------------------
plot(sim.ssn, "Sim_Values",
   xlab = "x-coordinate", ylab = "y-coordinate",
   cex = 1.5)

## ----eval=FALSE----------------------------------------------------------
## plot(sim.ssn, "Sim_Values",
##    xlab = "x-coordinate", ylab = "y-coordinate",
##    cex = 1.5)

## ------------------------------------------------------------------------
simDFobs <- getSSNdata.frame(sim.ssn, "Obs")
simDFpred <- getSSNdata.frame(sim.ssn, "preds")

## ------------------------------------------------------------------------
simpreds <- simDFpred[,"Sim_Values"]
simDFpred[,"Sim_Values"] <- NA
sim.ssn <- putSSNdata.frame(simDFpred, sim.ssn, "preds")

## ----cache=FALSE---------------------------------------------------------
glmssn.out <- glmssn(Sim_Values ~ X1 + X2 + F1, sim.ssn,
   CorModels = c("LinearSill.tailup", "Mariah.taildown",
   "Exponential.Euclid", "RE1", "RE2"),
   addfunccol = "addfunccol")

## ------------------------------------------------------------------------
summary(glmssn.out)

## ----SimTvP,include=FALSE,cache=FALSE------------------------------------
glmssn.pred <- predict(glmssn.out,"preds")
predDF <- getSSNdata.frame(glmssn.pred, "preds")
plot(simpreds, predDF[,"Sim_Values"], xlab = "True",
   ylab = "Predicted", pch = 19)

## ----eval=FALSE,cache=FALSE----------------------------------------------
## glmssn.pred <- predict(glmssn.out,"preds")
## predDF <- getSSNdata.frame(glmssn.pred, "preds")
## plot(simpreds, predDF[,"Sim_Values"], xlab = "True",
##    ylab = "Predicted", pch = 19)

Try the SSN package in your browser

Any scripts or data that you put into this service are public.

SSN documentation built on March 7, 2023, 5:30 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close