## ----setup, include=FALSE, cache=FALSE------------------------------
require(knitr)
opts_knit$set(aliases=c(h = 'fig.height', w = 'fig.width'),
              out.format='latex', use.highlight=TRUE)
opts_chunk$set(fig.path='kgraph/exPPA-', fig.align='center', fig.show='hold',
               prompt=FALSE, fig.width=6, fig.height=6,
               out.width='.65\\linewidth', size='scriptsize',
               cache=FALSE, message=FALSE,
               tidy = TRUE)
options(replace.assign=TRUE, width=70)
options(warn = -1)


## -------------------------------------------------------------------
require(spatstat)
data(japanesepines)
class(japanesepines)
str(japanesepines)
summary(japanesepines)


## -------------------------------------------------------------------
plot.ppp(japanesepines, main="Locations of Japanese pine trees", axes=T); grid()


## -------------------------------------------------------------------
data(redwoodfull)
data(cells)


## ----three-ppp, out.width='\\linewidth'-----------------------------
par(mfrow=c(2,2))
plot.ppp(japanesepines, main="Japanese Pines", axes=T) 
plot.ppp(redwoodfull, main="Redwoods", axes=T) 
plot.ppp(cells, main="Biological cells", axes=T) 
par(mfrow=c(1,1))


## -------------------------------------------------------------------
G <- Gest(japanesepines)
class(G)
summary(G)


## -------------------------------------------------------------------
plot(G, main="G-function, Japanese pines")


## -------------------------------------------------------------------
plot(G, xlim=c(0,.16), main="G-function, Japanese pines")


## -------------------------------------------------------------------
1 - exp(-japanesepines$n * pi * (0.16)^2)


## -------------------------------------------------------------------
G$r[min(which(G$rs >= 0.95))]
G$rs[min(which(G$r >= 0.05))]


## ----show-random, tidy=FALSE----------------------------------------
## 1 : example of random pattern
plot(G, cbind(rs,theo) ~ theo,
     main="G-function, Japanese pines") 


## ----G-redwoods, fig.width=10, fig.height=5, tidy=FALSE, out.width='\\linewidth'----
## 2 : example of clustered pattern
G.clust <- Gest(redwoodfull)
par(mfrow=c(1,2))
plot(G.clust, cbind(rs,theo) ~ r,
     main="G-function, Redwood trees")
plot(G.clust, cbind(rs,theo) ~ theo,
     main="G-function, Redwood trees")
par(mfrow=c(1,1))


## ----G-cells, fig.width=12, fig.height=6, out.width='\\linewidth', tidy=FALSE----
## 3 : example of dispersed / regular pattern
G.disp <- Gest(cells)
par(mfrow=c(1,2))
plot(G.disp, cbind(rs,theo) ~ r,
     main="G-function, cells")
plot(G.disp, cbind(rs,theo) ~ theo,
     main="G-function, cells")
par(mfrow=c(1,1))


## ----G-jap, tidy=FALSE----------------------------------------------
set.seed(30)
rmax.jap <- G$r[min(which(G$rs > 0.98))]
r <- seq(0, rmax.jap, by = 0.005)
envjap <- envelope(japanesepines, fun=Gest,
                   r=r, nrank=2, nsim=99, verbose=F)
plot(envjap, xlim=c(0, rmax.jap),
     main="Japanese pines, G-function envelope")


## ----G-env-1, fig.width=12, fig.height=6, out.width='\\linewidth', tidy=FALSE----
rmax.red <- G.clust$r[min(which(G.clust$rs > 0.98))]
r <- seq(0, rmax.red, by = 0.005)
envred <- envelope(redwoodfull, fun=Gest, r=r, nrank=2, nsim=99, verbose=F)
rmax.cell <- G.disp$r[min(which(G.disp$rs > 0.98))]
r <- seq(0, rmax.cell, by = 0.005)
envcells <- envelope(cells, fun=Gest, r=r, nrank=2, nsim=99, verbose=F)


## ----G-env-2, fig.width=12, fig.height=6, out.width='\\linewidth', tidy=FALSE----
par(mfrow=c(1,2))
plot(envred, xlim=c(0, rmax.red),
     main="Redwood trees, G-function envelope")
plot(envcells, xlim=c(0, rmax.cell),
     main="Cells, G-function envelope")
par(mfrow=c(1,1))


## ----optimum.bandwidth----------------------------------------------
print(bw.o <- bw.diggle(redwoodfull))
print(bw.CvL(redwoodfull))
print(bw.ppl(redwoodfull))


## ----compute-plot-densities, out.width='\\linewidth'----------------
d1 <- density.ppp(redwoodfull, sigma = bw.o, 
                  kernel = "quartic")
class(d1)
d05 <- density.ppp(redwoodfull, sigma = bw.o, 
                   kernel = "quartic", adjust = 0.5)
d4 <- density.ppp(redwoodfull, sigma = bw.o, 
                   kernel = "quartic", adjust = 4)
d2 <- density.ppp(redwoodfull, sigma = bw.o, 
                  kernel = "quartic", adjust = 2)
par(mfrow=c(2,2))
plot.im(d05, main=paste("Bandwidth=",round(bw.o*.4, 4),
                     " (optimum*.5)", sep="", collapse=""))
contour(d05, add = TRUE)
plot(d1, main=paste("Bandwidth=", round(bw.o, 4), 
                    " (optimum)", sep="", collapse=""))
contour(d1, add = TRUE)
plot(d2, main=paste("Bandwidth=",round(bw.o*1.5, 4), 
                     " (optimum*1.5)", sep="", collapse=""))
contour(d2, add = TRUE)
plot(d4, main=paste("Bandwidth=", round(bw.o*2,4),
                    " (optimum*2)", sep="", collapse=""))
contour(d4, add = TRUE)
par(mfrow=c(1,1))


## ----K-1, out.width='.8\\linewidth', tidy=FALSE---------------------
Kjap <- Kest(japanesepines)
Kred <- Kest(redwoodfull)
Kcells <- Kest(cells)


## ----K-2, fig.width=12, fig.height=4, out.width='\\linewidth', tidy=FALSE----
par(mfrow=c(1,3))
plot(Kjap, main="Japanese pines")
plot(Kred, main="Redwoods")
plot(Kcells, main="Cells")
par(mfrow=c(1,1))


## ----K-env, fig.width=12, fig.height=4, out.width='\\linewidth', tidy=FALSE----
r <- seq(0, sqrt(2)/6, by = 0.005)
envjap <- envelope(japanesepines, fun=Kest, r=r, nrank=2, nsim=99, verbose=F)
envred <- envelope(redwoodfull, fun=Kest, r=r, nrank=2, nsim=99, verbose=F)
envcells <- envelope(cells, fun=Kest, r=r, nrank=2, nsim=99, verbose=F)


## ----fig.width=12, fig.height=4, out.width='\\linewidth'------------
par(mfrow=c(1,3))
plot(envjap, main="Japanese pines, K-function envelope")
plot(envred, main="Redwood trees, K-function envelope")
plot(envcells, main="Cells, K-function envelope")
par(mfrow=c(1,1))


## -------------------------------------------------------------------
r <- seq(0, sqrt(2)/6, by = 0.005)
envjap <- envelope(japanesepines, fun=Lest, r=r, nrank=2, nsim=99, verbose=F)
envred <- envelope(redwoodfull, fun=Lest, r=r, nrank=2, nsim=99, verbose=F)
envcells <- envelope(cells, fun=Lest, r=r, nrank=2, nsim=99, verbose=F)


## ----fig.width=12, fig.height=4, out.width='\\linewidth'------------
par(mfrow=c(1,3))
plot(envjap, main="Japanese pines, L-function envelope")
plot(envred, main="Redwood trees, L-function envelope")
plot(envcells, main="Cells, L-function envelope")
par(mfrow=c(1,1))


## ----subset-windows-------------------------------------------------
str(japanesepines, max.level=1)
print(japanesepines$window)
(window.nw <- owin(xrange=c(0,0.5), yrange=c(0.5,1)))
table(is.in <- inside.owin(japanesepines,  w=window.nw))
japanesepines.nw <- japanesepines[is.in]


## ----show-window----------------------------------------------------
Window(japanesepines.nw)


## ----reduce-window--------------------------------------------------
Window(japanesepines.nw) <- window.nw
Window(japanesepines.nw)


## ----fig.width=5, fig.height=5--------------------------------------
plot(japanesepines.nw, main="Japanese pines, NW quadrant")


## ----fig.width=12, fig.height=12, out.width='\\linewidth'-----------
par(mfrow=c(2,2))
G <- Gest(japanesepines.nw)
plot(G, main = "G-function, Japanese pines, NW quadrant", xlim=c(0,.1))
G <- Gest(japanesepines)
plot(G, main = "G-function, Japanese pines, all", xlim=c(0,.1))
L <- Lest(japanesepines.nw)
plot(L, main = "L-function, Japanese pines, NW quadrant", xlim=c(0,.1))
L <- Lest(japanesepines)
plot(L, main = "L-function, Japanese pines, all", xlim=c(0,.1))
par(mfrow=c(1,1))


## ----tidy=FALSE, fig.width=10, fig.height=10, out.width='0.8\\linewidth'----
data(meuse, package = "sp")
meuse <- meuse[,c("x","y","ffreq")]
require(sf)
meuse.sf <- st_as_sf(meuse, coords = 1:2)
meuse.ppp <- as.ppp(meuse.sf)
str(meuse.ppp)
marks(meuse.ppp) <- meuse$ffreq
tmp <- plot(meuse.ppp, use.marks=TRUE, 
     cols=c("red","orange","green"), 
     chars=16, which.marks="marks", 
     main="Meuse floodplain flood frequency class", 
     axes=T)
grid()
legend("left", pch=16, 
       col=c("red","orange","green"), 
       legend=c("Annually","2-5 Years", "> 5 Years"))


## -------------------------------------------------------------------
meuse.ppp$window


## ----tidy=FALSE, fig.width=10, fig.height=10, out.width='0.8\\linewidth'----
meuse.ppp.r <- meuse.ppp
(meuse.ppp.r$window <- ripras(meuse.ppp))
tmp <- plot(meuse.ppp.r, use.marks=TRUE, 
     cols=c("red","orange","green"), 
     chars=16, which.marks="ffreq",
     main="Meuse floodplain flood frequency class", 
     boundary=2, axes=T)
grid()
legend("left", pch=16, 
       col=c("red","orange","green"), 
       legend=c("Annually","2-5 Years", "> 5 Years"))
ch <- convexhull(meuse.ppp)
lines(ch$bdry[[1]]$x, ch$bdry[[1]]$y, lty=2)
legend("bottomright", lty=1:2, 
       legend=c("Ripley-Rasson", "convex hull"))


## -------------------------------------------------------------------
str(meuse.ppp)
intensity(meuse.ppp)
intensity(meuse.ppp.r)
round(intensity(meuse.ppp.r)/intensity(meuse.ppp),2)


## ----fig.width=12, fig.height=6, out.width='\\linewidth'------------
par(mfrow=c(1,2))
plot(Fest(meuse.ppp), main="rectangular window", xlim=c(0,550))
plot(Fest(meuse.ppp.r), main="polygonal window", xlim=c(0,550))
par(mfrow=c(1,2))


## ----fig.width=8, fig.height=8, out.width='\\linewidth'-------------
data(clmfires)
plot(clmfires, which.marks="cause", bg=2:5, chars=21:24, cex=0.5, axes=T,
     main="Castilla-La Mancha forest fires")
grid()
legend("topleft", pch=21:24, legend=levels(clmfires$marks$cause), pt.bg=2:5)


## ----climfires-split, fig.width=10, fig.height=10, out.width='\\linewidth', tidy=FALSE----
clmfires.split <- split(clmfires)
str(clmfires.split, max.level=1)
plot(clmfires.split, use.marks=FALSE,
     main="Castilla-La Mancha forest fires",
     pch=21, bg=2)


## -------------------------------------------------------------------
str(clmfires$marks)


## -------------------------------------------------------------------
clmfires.cause <- clmfires
is.multitype(clmfires.cause)
clmfires.cause$marks <- clmfires$marks$cause
is.multitype(clmfires.cause)


## -------------------------------------------------------------------
Kcross.il <- Kcross(clmfires.cause, "intentional", "lightning", correction="translate")
plot(Kcross.il); grid()


## ----superimpose, out.width='0.8\\linewidth'------------------------
two.trees <- superimpose(rw=redwoodfull, jp=japanesepines)
str(two.trees)
plot(two.trees, main="Superimposed point patterns", cols=c("green","blue"))


## ----kcross-jp-rw---------------------------------------------------
Kcross.jp.rw <- Kcross(two.trees)
plot(Kcross.jp.rw); grid()


## ----show-longleaf, out.width='0.8\\linewidth'----------------------
data(longleaf)
summary(longleaf)
Window(longleaf)
plot(longleaf, main="Longleaf pines, location and DBH",
     cols=function(x) ifelse(x < 30, "green", "red"))
grid()


## ----longleaf-K-----------------------------------------------------
K.long <- Kest(longleaf)
plot(K.long, main="K function, longleaf pines")


## ----out.width='0.8\\linewidth'-------------------------------------
clmfires.i <- split(clmfires, "cause")$intentional
plot(clmfires.i, chars=21, cex=0.5, bg=2, axes=T, main="Castilla-La Mancha intentional forest fires", use.marks=FALSE)
grid()


## -------------------------------------------------------------------
marks(clmfires.i) <- NULL


## -------------------------------------------------------------------
print(m.pois <- ppm(clmfires.i, trend=~1, interaction=NULL))
class(m.pois)
exp(coef(m.pois))
intensity(clmfires.i)
(clmfires.i$n/summary(clmfires.i)$window$areas)


## -------------------------------------------------------------------
(m.ts1 <- ppm(clmfires.i, trend=~polynom(x,y,1), interaction=NULL))
(m.ts2 <- ppm(clmfires.i, trend=~polynom(x,y,2), interaction=NULL))


## ----fig.width=12, fig.height=6, out.width='\\linewidth'------------
par(mfrow=c(1,2))
plot.ppm(m.ts1, ngrid=c(80,80), how="image", superimpose=F, trend=T, se=F, pause=F, main="1st-order trend")
plot.ppm(m.ts2, ngrid=c(80,80), how="image", superimpose=F, trend=T, se=F, pause=F, main="2nd-order trend")
par(mfrow=c(1,1))


## ----fig.width=12, fig.height=6, out.width='\\linewidth'------------
par(mfrow=c(1,2))
plot.ppm(m.ts1, ngrid=c(80,80), how="persp", theta=-30, phi=30, trend=T, se=F, pause=F, main="1st-order trend")
plot.ppm(m.ts2, ngrid=c(80,80), how="persp", theta=-30, phi=30, trend=T, se=F, pause=F, main="2nd-order trend")
par(mfrow=c(1,1))


## -------------------------------------------------------------------
anova.ppm(m.ts2, m.ts1, m.pois)


## ----plot-Kest-climfires-i------------------------------------------
plot(Kest(clmfires.i, r=seq(0,10,by=.2)))


## ----strauss.4km, tidy=FALSE----------------------------------------
(m.strauss.4 <- ppm(clmfires.i, trend=~1,
                    interaction=StraussHard(r=4)))
exp(coef(m.strauss.4))


## ----compare-logLik, tidy=FALSE-------------------------------------
data.frame(
    model=c("Poisson","1st order trend", "2nd order trend", "Strauss/hard core"),
    likelihood=c(logLik(m.pois, warn=F),
        logLik(m.ts1, warn=F),
        logLik(m.ts2, warn=F),
        logLik(m.strauss.4, warn=F)))


## ----show.climfires-------------------------------------------------
names(clmfires.extra)
names(clmfires.extra$clmcov200)
names(clmfires.extra$clmcov200$landuse)
names(clmfires.extra$clmcov200$landuse)
levels(clmfires.extra$clmcov200$landuse$v)


## ----out.width='\\linewidth'----------------------------------------
plot(clmfires.extra$clmcov200, main="200 m grid covariates")


## ----fig.align='center'---------------------------------------------
require(terra)
clmfires.lu.grid <- rast(clmfires.extra$clmcov200$landuse)
class(clmfires.lu.grid)
plot(clmfires.lu.grid, col=hcl.colors(11, palette = "viridis"))


## ----linear.in.covars, tidy=FALSE-----------------------------------
levels(clmfires.extra$clmcov200$landuse)
(m.lu <- ppm(clmfires.i, ~ landuse - 1,
             covariates=clmfires.extra$clmcov200,
             interaction=NULL))
data.frame(model=c("Poisson", "Landuse"),
           likelihood=c(logLik(m.pois),
               logLik(m.lu, warn=F)))


## ----strauss.in.covars, tidy=FALSE----------------------------------
(m.lu.strauss.4 <- ppm(clmfires.i, ~ landuse,
                       covariates=clmfires.extra$clmcov200,
                       interaction=StraussHard(r = 4)))
data.frame(
    model=c("Poisson", "Landuse", "Strauss", "Landuse + Strauss"),
    likelihood=c(logLik(m.pois),
        logLik(m.lu, warn=F),
        logLik(m.strauss.4, warn=F),
        logLik(m.lu.strauss.4, warn=F)))


## ----predict-ppm-lu, out.width='0.8\\linewidth'---------------------
pred.lu <- predict(m.lu, covariates=clmfires.extra$clmcov200)
summary(pred.lu)
image(pred.lu, main="Fire intensity based on land use")


## ----predict-ppm-lu-strauss, tidy=FALSE-----------------------------
pred.lu.strauss <- predict(m.lu.strauss.4,
                           covariates=clmfires.extra$clmcov200)
summary(pred.lu.strauss)
range(pred.lu.strauss)
range(pred.lu)
mean(pred.lu.strauss)
mean(pred.lu)


## ----plot-predict-ppm-lu-strauss, tidy=FALSE, fig.width=12, fig.height=12, out.width='\\linewidth'----
par(mfrow=c(1,2))
zlim=c(min(pred.lu, pred.lu.strauss), 
       max(pred.lu, pred.lu.strauss))
image(pred.lu, zlim=zlim,
      main="land use")
image(pred.lu.strauss, zlim=zlim,
      main="land use + Strauss process")
par(mfrow=c(1,2))


## ----load-stpp------------------------------------------------------
library("stpp")
data("fmd")
data("northcumbria")
summary(fmd)
dim(fmd)


## ----help-stpp, eval=FALSE------------------------------------------
## help(fmd)


## ----hist-stpp, tidy=FALSE------------------------------------------
hist(fmd[,"ReportedDay"], xlab="reported day", 
     main="Cases of Foot-and-mouth disease", breaks=16)
rug(fmd[,"ReportedDay"])


## ----convert-stpp---------------------------------------------------
class(fmd)
fmd <- as.3dpoints(fmd)
class(fmd)


## ----plot-stpp-locations, out.width = '\\linewidth'-----------------
plot(fmd, s.region = northcumbria, pch = 21, mark = TRUE,
     mark.col=0, mark.cexmin=0.2, mark.cexmax=1.2, col="blue", bg="red")


## ----slice-stpp-----------------------------------------------------
dim(fmd)[1]
fmd.1 <- as.3dpoints(fmd[fmd[,3] <= 50,])
fmd.2 <- as.3dpoints(fmd[(fmd[,3] > 50) & (fmd[,3] <= 100),])
fmd.3 <- as.3dpoints(fmd[(fmd[,3] > 100) & (fmd[,3] <= 150),])
fmd.4 <- as.3dpoints(fmd[fmd[,3] > 150,])
dim(fmd.1)[1]; dim(fmd.2)[1]; dim(fmd.3)[1]; dim(fmd.4)[1]


## ----ncumbria-slice1, out.width = '0.8\\linewidth'------------------
plot(fmd.1, s.region = northcumbria, pch=21, col="blue", bg="red", mark=T, mark.col=0, mark.cexmin=1, mark.cexmax=1)
title("Days 0-50")
grid()


## ----ncumbria-slice2, out.width = '0.8\\linewidth'------------------
plot(fmd.2, s.region = northcumbria, pch=21, col="blue", bg="red", mark=T, mark.col=0, mark.cexmin=1, mark.cexmax=1)
title("Days 51-100")
grid()


## ----ncumbria-slice3, out.width = '0.8\\linewidth'------------------
plot(fmd.3, s.region = northcumbria, pch=21, col="blue", bg="red", mark=T, mark.col=0, mark.cexmin=1, mark.cexmax=1)
title("Days 101-150")
grid()


## ----ncumbria-slice4, out.width = '0.8\\linewidth'------------------
plot(fmd.4, s.region = northcumbria, pch=21, col="blue", bg="red", mark=T, mark.col=0, mark.cexmin=1, mark.cexmax=1)
title("Days 151-200")
grid()


## ----compute-owin-stpp----------------------------------------------
w <- owin(poly=list(x=northcumbria[,1], y=northcumbria[,2]))


## ----intensity-stpp-------------------------------------------------
1/(intensity(fmd.1.ppp <- ppp(fmd.1[,1], fmd.1[,2], window=w))*10^6)
1/(intensity(fmd.2.ppp <- ppp(fmd.2[,1], fmd.2[,2], window=w))*10^6)
1/(intensity(fmd.3.ppp <- ppp(fmd.3[,1], fmd.3[,2], window=w))*10^6)
1/(intensity(fmd.4.ppp <- ppp(fmd.4[,1], fmd.4[,2], window=w))*10^6)


## ----Fest-stpp, fig.width=12, fig.height=12, out.width='\\linewidth'----
par(mfrow=c(2,2))
plot(Fest(fmd.1.ppp), main="Days 0-50")
plot(Fest(fmd.2.ppp), main="Days 51-100")
plot(Fest(fmd.3.ppp), main="Days 101-150")
plot(Fest(fmd.4.ppp), main="Days 151-200")
par(mfrow=c(1,1))


## ----Gest-stpp, fig.width=12, fig.height=12, out.width='\\linewidth'----
par(mfrow=c(2,2))
plot(Gest(fmd.1.ppp), main="Days 0-50")
plot(Gest(fmd.2.ppp), main="Days 51-100")
plot(Gest(fmd.3.ppp), main="Days 101-150")
plot(Gest(fmd.4.ppp), main="Days 151-200")
par(mfrow=c(1,1))


## ----marked-stpp, out.width = '\\linewidth'-------------------------
fmd.all.ppp <- superimpose(Q1=fmd.1.ppp, Q2=fmd.2.ppp, Q3=fmd.3.ppp, Q4=fmd.4.ppp)
plot(fmd.all.ppp, main="2001 Foot-and-mouth disease, 50-day intervals", 
     cex=0.9, pch=21, col=1, bg=2:5)


## ----kcross-stpp-compute--------------------------------------------
Kcross.1.2 <- Kcross(fmd.all.ppp, "Q1", "Q2")
Kcross.2.3 <- Kcross(fmd.all.ppp, "Q2", "Q3")
Kcross.3.4 <- Kcross(fmd.all.ppp, "Q3", "Q4")


## ----kcross-stpp, fig.width=12, fig.height=4, out.width='\\linewidth'----
par(mfrow=c(1,3))
plot(Kcross.1.2, main="0-50 vs.\ 51-100")
plot(Kcross.2.3, main="51-100 vs.\ 101-150")
plot(Kcross.3.4, main="101-150 vs.\ 151-200")
par(mfrow=c(1,1))


## ----readBaltimore--------------------------------------------------
baltim <- sf::st_read(system.file("shapes/baltim.shp", 
                               package="maptools"))
str(baltim)
st_crs(baltim)  # no coordinate system is defined for this shapefile


## ----convertBaltimore, tidy=FALSE-----------------------------------
require(spatstat)
baltim.ppp <- as.ppp(baltim)
str(baltim.ppp)
window(baltim.ppp)
marks(baltim.ppp) <- st_drop_geometry(baltim)
summary(marks(baltim.ppp))
#
op <- par(no.readonly = TRUE)
par(mar=rep(0.5, 4))
plot(baltim.ppp, which.marks = "PRICE", axes=T, 
     main="Baltimore house sale prices, k$")
par(op)


## ----exportMeuse----------------------------------------------------
data(meuse, package = "sp")
class(meuse)
write.csv(meuse, file="tmp.csv", row.names = FALSE)


## ----eval=FALSE-----------------------------------------------------
## file.show("tmp.csv")


## ----importMeuse----------------------------------------------------
pp <- read.csv(file="tmp.csv")
class(pp)
pp$ffreq <- as.factor(pp$ffreq)
pp$soil <- as.factor(pp$soil)
pp$lime <- as.factor(pp$lime)
str(pp)


## ----convertMeuseSf-------------------------------------------------
require(sf)
pp.sf <- st_as_sf(pp, coords = c("x", "y"), crs = 28992)
class(pp.sf)
summary(pp.sf)


## ----convertMeusePpp------------------------------------------------
pp.m <- as.ppp(pp.sf)
class(pp.m)
summary(pp.m)
window(pp.m)


## ----changeMarksMeusePpp--------------------------------------------
summary(marks(pp.m))
marks(pp.m) <- pp.sf$soil
table(marks(pp.m))


## ----showMeusePpp, tidy=FALSE---------------------------------------
op <- par(no.readonly = TRUE)
par(mar=rep(0.5, 4))
plot(pp.m,  cols=c("red","blue","green"),
     pch=20, main="Meuse soil types", axes=T)
par(op)


## ----read-meuse-boundary-not, eval=FALSE----------------------------
## meuseBoundary <- sf::st_read(dsn=".", layer="meuseBoundary")
## class(meuseBoundary)


## ----read-meuse-boundary, echo=FALSE--------------------------------
meuseBoundary <- sf::st_read("./ds/meuse", "meuseBoundary")
class(meuseBoundary)


## ----show-boundary--------------------------------------------------
class(meuseBoundary)
(meuseBoundary.win <- as.owin(meuseBoundary))


## ----set-boundary-window--------------------------------------------
pp.m$window <- meuseBoundary.win


## ----plotMeusePppWithBoundary, tidy=FALSE---------------------------
op <- par(no.readonly = TRUE)
par(mar=rep(0.5, 4))
plot(pp.m, cols=c("red","blue","green"),
     pch=20, main="Meuse soil types", axes=T)
par(op)


## -------------------------------------------------------------------
unlink("tmp.csv")
rm(meuse, pp, pp.sf, pp.m, meuseBoundary, meuseBoundary.win, op)